neuper@37906: (* rationals, i.e. fractions of multivariate polynomials over the real field neuper@37906: author: isac team neuper@37906: Copyright (c) isac team 2002 neuper@37906: Use is subject to license terms. neuper@37906: neuper@37906: depends on Poly (and not on Atools), because neuper@37906: fractions with _normalized_ polynomials are canceled, added, etc. neuper@37906: *) neuper@37906: neuper@37950: theory Rational imports Poly begin neuper@37906: neuper@37906: consts neuper@37906: neuper@37906: is'_expanded :: "real => bool" ("_ is'_expanded") (*RL->Poly.thy*) neuper@37906: is'_ratpolyexp :: "real => bool" ("_ is'_ratpolyexp") neuper@37906: neuper@37950: axioms (*.not contained in Isabelle2002, neuper@37950: stated as axioms, TODO?: prove as theorems*) neuper@37950: neuper@37978: mult_cross: "[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)" neuper@37978: mult_cross1: " b ~= 0 ==> (a / b = c ) = (a = b * c)" neuper@37978: mult_cross2: " d ~= 0 ==> (a = c / d) = (a * d = c)" neuper@37950: neuper@37978: add_minus: "a + b - b = a"(*RL->Poly.thy*) neuper@37978: add_minus1: "a - b + b = a"(*RL->Poly.thy*) neuper@37950: neuper@37978: rat_mult: "a / b * (c / d) = a * c / (b * d)"(*?Isa02*) neuper@37978: rat_mult2: "a / b * c = a * c / b "(*?Isa02*) neuper@37978: neuper@37978: rat_mult_poly_l: "c is_polyexp ==> c * (a / b) = c * a / b" neuper@37978: rat_mult_poly_r: "c is_polyexp ==> (a / b) * c = a * c / b" neuper@37906: neuper@37906: (*real_times_divide1_eq .. Isa02*) neuper@37978: real_times_divide_1_eq: "-1 * (c / d) =-1 * c / d " neuper@37978: real_times_divide_num: "a is_const ==> neuper@37950: a * (c / d) = a * c / d " neuper@37906: neuper@37978: real_mult_div_cancel2: "k ~= 0 ==> m * k / (n * k) = m / n" neuper@37978: (*real_mult_div_cancel1: "k ~= 0 ==> k * m / (k * n) = m / n"..Isa02*) neuper@37906: neuper@37978: real_divide_divide1: "y ~= 0 ==> (u / v) / (y / z) = (u / v) * (z / y)" neuper@37978: real_divide_divide1_mg: "y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)" neuper@37978: (*real_divide_divide2_eq: "x / y / z = x / (y * z)"..Isa02*) neuper@37906: neuper@37978: rat_power: "(a / b)^^^n = (a^^^n) / (b^^^n)" neuper@37978: neuper@37978: neuper@37978: rat_add: "[| a is_const; b is_const; c is_const; d is_const |] ==> neuper@37950: a / c + b / d = (a * d + b * c) / (c * d)" neuper@37978: rat_add_assoc: "[| a is_const; b is_const; c is_const; d is_const |] ==> neuper@37950: a / c +(b / d + e) = (a * d + b * c)/(d * c) + e" neuper@37978: rat_add1: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a / c + b / c = (a + b) / c" neuper@37978: rat_add1_assoc: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a / c + (b / c + e) = (a + b) / c + e" neuper@37978: rat_add2: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a / c + b = (a + b * c) / c" neuper@37978: rat_add2_assoc: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a / c + (b + e) = (a + b * c) / c + e" neuper@37978: rat_add3: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a + b / c = (a * c + b) / c" neuper@37978: rat_add3_assoc: "[| a is_const; b is_const; c is_const |] ==> neuper@37950: a + (b / c + e) = (a * c + b) / c + e" neuper@37950: neuper@37950: text {*calculate in rationals: gcd, lcm, etc. neuper@37950: (c) Stefan Karnel 2002 neuper@37950: Institute for Mathematics D and Institute for Software Technology, neuper@37950: TU-Graz SS 2002 *} neuper@37950: neuper@37950: text {* Remark on notions in the documentation below: neuper@37950: referring to the remark on 'polynomials' in Poly.sml we use neuper@37950: [2] 'polynomial' normalform (Polynom) neuper@37950: [3] 'expanded_term' normalform (Ausmultiplizierter Term), neuper@37950: where normalform [2] is a special case of [3], i.e. [3] implies [2]. neuper@37950: Instead of neuper@37950: 'fraction with numerator and nominator both in normalform [2]' neuper@37950: 'fraction with numerator and nominator both in normalform [3]' neuper@37950: we say: neuper@37950: 'fraction in normalform [2]' neuper@37950: 'fraction in normalform [3]' neuper@37950: or neuper@37950: 'fraction [2]' neuper@37950: 'fraction [3]'. neuper@37950: a 'simple fraction' is a term with '/' as outmost operator and neuper@37950: numerator and nominator in normalform [2] or [3]. neuper@37950: *} neuper@37950: neuper@37979: ML {* neuper@37972: val thy = @{theory}; neuper@37972: neuper@37950: signature RATIONALI = neuper@37950: sig neuper@37950: type mv_monom neuper@37950: type mv_poly neuper@37950: val add_fraction_ : theory -> term -> (term * term list) option neuper@37950: val add_fraction_p_ : theory -> term -> (term * term list) option neuper@37950: val calculate_Rational : rls neuper@37950: val calc_rat_erls:rls neuper@37950: val cancel : rls neuper@37950: val cancel_ : theory -> term -> (term * term list) option neuper@37950: val cancel_p : rls neuper@37950: val cancel_p_ : theory -> term -> (term * term list) option neuper@37950: val common_nominator : rls neuper@37950: val common_nominator_ : theory -> term -> (term * term list) option neuper@37950: val common_nominator_p : rls neuper@37950: val common_nominator_p_ : theory -> term -> (term * term list) option neuper@37950: val eval_is_expanded : string -> 'a -> term -> theory -> neuper@37950: (string * term) option neuper@37950: val expanded2polynomial : term -> term option neuper@37950: val factout_ : theory -> term -> (term * term list) option neuper@37950: val factout_p_ : theory -> term -> (term * term list) option neuper@37950: val is_expanded : term -> bool neuper@37950: val is_polynomial : term -> bool neuper@37950: neuper@37950: val mv_gcd : (int * int list) list -> mv_poly -> mv_poly neuper@37950: val mv_lcm : mv_poly -> mv_poly -> mv_poly neuper@37950: neuper@37950: val norm_expanded_rat_ : theory -> term -> (term * term list) option neuper@37950: (*WN0602.2.6.pull into struct !!! neuper@37950: val norm_Rational : rls(*.normalizes an arbitrary rational term without neuper@37950: roots into a simple and canceled fraction neuper@37950: with normalform [2].*) neuper@37950: *) neuper@37950: (*val norm_rational_p : 19.10.02 missing FIXXXXXXXXXXXXME neuper@37950: rls (*.normalizes an rational term [2] without neuper@37950: roots into a simple and canceled fraction neuper@37950: with normalform [2].*) neuper@37950: *) neuper@37950: val norm_rational_ : theory -> term -> (term * term list) option neuper@37950: val polynomial2expanded : term -> term option neuper@37950: val rational_erls : neuper@37950: rls (*.evaluates an arbitrary rational term with numerals.*) neuper@37950: neuper@37950: (*WN0210???SK: fehlen Funktionen, die exportiert werden sollen ? *) neuper@37906: end neuper@37950: neuper@37950: (*.************************************************************************** neuper@37950: survey on the functions neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~ neuper@37950: [2] 'polynomial' :rls | [3]'expanded_term':rls neuper@37950: --------------------:------------------+-------------------:----------------- neuper@37950: factout_p_ : | factout_ : neuper@37950: cancel_p_ : | cancel_ : neuper@37950: :cancel_p | :cancel neuper@37950: --------------------:------------------+-------------------:----------------- neuper@37950: common_nominator_p_: | common_nominator_ : neuper@37950: :common_nominator_p| :common_nominator neuper@37950: add_fraction_p_ : | add_fraction_ : neuper@37950: --------------------:------------------+-------------------:----------------- neuper@37950: ???SK :norm_rational_p | :norm_rational neuper@37950: neuper@37950: This survey shows only the principal functions for reuse, and the identifiers neuper@37950: of the rls exported. The list below shows some more useful functions. neuper@37950: neuper@37950: neuper@37950: conversion from Isabelle-term to internal representation neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ neuper@37950: neuper@37950: ... BITTE FORTSETZEN ... neuper@37950: neuper@37950: polynomial2expanded = ... neuper@37950: expanded2polynomial = ... neuper@37950: neuper@37950: remark: polynomial2expanded o expanded2polynomial = I, neuper@37950: where 'o' is function chaining, and 'I' is identity WN0210???SK neuper@37950: neuper@37950: functions for greatest common divisor and canceling neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ neuper@37950: mv_gcd neuper@37950: factout_ neuper@37950: factout_p_ neuper@37950: cancel_ neuper@37950: cancel_p_ neuper@37950: neuper@37950: functions for least common multiple and addition of fractions neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ neuper@37950: mv_lcm neuper@37950: common_nominator_ neuper@37950: common_nominator_p_ neuper@37950: add_fraction_ (*.add 2 or more fractions.*) neuper@37950: add_fraction_p_ (*.add 2 or more fractions.*) neuper@37950: neuper@37950: functions for normalform of rationals neuper@37950: ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ neuper@37950: WN0210???SK interne Funktionen f"ur norm_rational: neuper@37950: schaffen diese SML-Funktionen wirklich ganz allgemeine Terme ? neuper@37950: neuper@37950: norm_rational_ neuper@37950: norm_expanded_rat_ neuper@37950: neuper@37950: **************************************************************************.*) neuper@37950: neuper@37950: neuper@37950: (*##*) neuper@37950: structure RationalI : RATIONALI = neuper@37950: struct neuper@37950: (*##*) neuper@37950: neuper@37950: infix mem ins union; (*WN100819 updating to Isabelle2009-2*) neuper@37950: fun x mem [] = false neuper@37950: | x mem (y :: ys) = x = y orelse x mem ys; neuper@37950: fun (x ins xs) = if x mem xs then xs else x :: xs; neuper@37950: fun xs union [] = xs neuper@37950: | [] union ys = ys neuper@37950: | (x :: xs) union ys = xs union (x ins ys); neuper@37950: neuper@37950: (*. gcd of integers .*) neuper@37950: (* die gcd Funktion von Isabelle funktioniert nicht richtig !!! *) neuper@37950: fun gcd_int a b = if b=0 then a neuper@37950: else gcd_int b (a mod b); neuper@37950: neuper@37950: (*. univariate polynomials (uv) .*) neuper@37950: (*. univariate polynomials are represented as a list neuper@37950: of the coefficent in reverse maximum degree order .*) neuper@37950: (*. 5 * x^5 + 4 * x^3 + 2 * x^2 + x + 19 => [19,1,2,4,0,5] .*) neuper@37950: type uv_poly = int list; neuper@37950: neuper@37950: (*. adds two uv polynomials .*) neuper@37950: fun uv_mod_add_poly ([]:uv_poly,p2:uv_poly) = p2:uv_poly neuper@37950: | uv_mod_add_poly (p1,[]) = p1 neuper@37950: | uv_mod_add_poly (x::p1,y::p2) = (x+y)::(uv_mod_add_poly(p1,p2)); neuper@37950: neuper@37950: (*. multiplies a uv polynomial with a skalar s .*) neuper@37950: fun uv_mod_smul_poly ([]:uv_poly,s:int) = []:uv_poly neuper@37950: | uv_mod_smul_poly (x::p,s) = (x*s)::(uv_mod_smul_poly(p,s)); neuper@37950: neuper@37950: (*. calculates the remainder of a polynomial divided by a skalar s .*) neuper@37950: fun uv_mod_rem_poly ([]:uv_poly,s) = []:uv_poly neuper@37950: | uv_mod_rem_poly (x::p,s) = (x mod s)::(uv_mod_smul_poly(p,s)); neuper@37950: neuper@37950: (*. calculates the degree of a uv polynomial .*) neuper@37950: fun uv_mod_deg ([]:uv_poly) = 0 neuper@37950: | uv_mod_deg p = length(p)-1; neuper@37950: neuper@37950: (*. calculates the remainder of x/p and represents it as neuper@37950: value between -p/2 and p/2 .*) neuper@37950: fun uv_mod_mod2(x,p)= neuper@37950: let neuper@37950: val y=(x mod p); neuper@37950: in neuper@37950: if (y)>(p div 2) then (y)-p else neuper@37950: ( neuper@37950: if (y)<(~p div 2) then p+(y) else (y) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*.calculates the remainder for each element of a integer list divided by p.*) neuper@37950: fun uv_mod_list_modp [] p = [] neuper@37950: | uv_mod_list_modp (x::xs) p = (uv_mod_mod2(x,p))::(uv_mod_list_modp xs p); neuper@37950: neuper@37950: (*. appends an integer at the end of a integer list .*) neuper@37950: fun uv_mod_null (p1:int list,0) = p1 neuper@37950: | uv_mod_null (p1:int list,n1:int) = uv_mod_null(p1,n1-1) @ [0]; neuper@37950: neuper@37950: (*. uv polynomial division, result is (quotient, remainder) .*) neuper@37950: (*. only for uv_mod_divides .*) neuper@37950: (* FIXME: Division von x^9+x^5+1 durch x-1000 funktioniert nicht, neuper@37950: integer zu klein *) neuper@37950: fun uv_mod_pdiv (p1:uv_poly) ([]:uv_poly) = neuper@37950: raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero") neuper@37950: | uv_mod_pdiv p1 [x] = neuper@37950: let neuper@38006: val xs= Unsynchronized.ref []; neuper@37950: in neuper@37950: if x<>0 then neuper@37950: ( neuper@37950: xs:=(uv_mod_rem_poly(p1,x)); neuper@37950: while length(!xs)>0 andalso hd(!xs)=0 do xs:=tl(!xs) neuper@37950: ) neuper@37950: else raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: division by zero"); neuper@37950: ([]:uv_poly,!xs:uv_poly) neuper@37950: end neuper@37950: | uv_mod_pdiv p1 p2 = neuper@37950: let neuper@37950: val n= uv_mod_deg(p2); neuper@38006: val m= Unsynchronized.ref (uv_mod_deg(p1)); neuper@38006: val p1'= Unsynchronized.ref (rev(p1)); neuper@37950: val p2'=(rev(p2)); neuper@37950: val lc2=hd(p2'); neuper@38006: val q= Unsynchronized.ref []; neuper@38006: val c= Unsynchronized.ref 0; neuper@38006: val output= Unsynchronized.ref ([],[]); neuper@37950: in neuper@37950: ( neuper@37950: if (!m)=0 orelse p2=[0] neuper@37950: then raise error ("RATIONALS_UV_MOD_PDIV_EXCEPTION: Division by zero") neuper@37950: else neuper@37950: ( neuper@37950: if (!m)=n do neuper@37950: ( neuper@37950: c:=hd(!p1') div hd(p2'); neuper@37950: if !c<>0 then neuper@37950: ( neuper@37950: p1':=uv_mod_add_poly(!p1',uv_mod_null(uv_mod_smul_poly(p2',~(!c)),!m-n)); neuper@37950: while length(!p1')>0 andalso hd(!p1')=0 do p1':= tl(!p1'); neuper@37950: m:=uv_mod_deg(!p1') neuper@37950: ) neuper@37950: else m:=0 neuper@37950: ); neuper@37950: output:=(rev(!q),rev(!p1')) neuper@37950: ) neuper@37950: ); neuper@37950: !output neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. divides p1 by p2 in Zp .*) neuper@37950: fun uv_mod_pdivp (p1:uv_poly) (p2:uv_poly) p = neuper@37950: let neuper@37950: val n=uv_mod_deg(p2); neuper@38006: val m= Unsynchronized.ref (uv_mod_deg(uv_mod_list_modp p1 p)); neuper@38006: val p1'= Unsynchronized.ref (rev(p1)); neuper@37950: val p2'=(rev(uv_mod_list_modp p2 p)); neuper@37950: val lc2=hd(p2'); neuper@38006: val q= Unsynchronized.ref []; neuper@38006: val c= Unsynchronized.ref 0; neuper@38006: val output= Unsynchronized.ref ([],[]); neuper@37950: in neuper@37950: ( neuper@37950: if (!m)=0 orelse p2=[0] then raise error ("RATIONALS_UV_MOD_PDIVP_EXCEPTION: Division by zero") neuper@37950: else neuper@37950: ( neuper@37950: if (!m)=n do neuper@37950: ( neuper@37950: c:=uv_mod_mod2(hd(!p1')*(power lc2 1), p); neuper@37950: q:=(!c)::(!q); neuper@37950: p1':=uv_mod_list_modp(tl(uv_mod_add_poly(uv_mod_smul_poly(!p1',lc2), neuper@37950: uv_mod_smul_poly(uv_mod_smul_poly(p2',hd(!p1')),~1)))) p; neuper@37950: m:=(!m)-1 neuper@37950: ); neuper@37950: neuper@37950: while !p1'<>[] andalso hd(!p1')=0 do neuper@37950: ( neuper@37950: p1':=tl(!p1') neuper@37950: ); neuper@37950: neuper@37950: output:=(rev(uv_mod_list_modp (!q) (p)),rev(!p1')) neuper@37950: ) neuper@37950: ); neuper@37950: !output:uv_poly * uv_poly neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the remainder of p1/p2 .*) neuper@37950: fun uv_mod_prest (p1:uv_poly) ([]:uv_poly) = raise error("UV_MOD_PREST_EXCEPTION: Division by zero") neuper@37950: | uv_mod_prest [] p2 = []:uv_poly neuper@37950: | uv_mod_prest p1 p2 = (#2(uv_mod_pdiv p1 p2)); neuper@37950: neuper@37950: (*. calculates the remainder of p1/p2 in Zp .*) neuper@37950: fun uv_mod_prestp (p1:uv_poly) ([]:uv_poly) p= raise error("UV_MOD_PRESTP_EXCEPTION: Division by zero") neuper@37950: | uv_mod_prestp [] p2 p= []:uv_poly neuper@37950: | uv_mod_prestp p1 p2 p = #2(uv_mod_pdivp p1 p2 p); neuper@37950: neuper@37950: (*. calculates the content of a uv polynomial .*) neuper@37950: fun uv_mod_cont ([]:uv_poly) = 0 neuper@37950: | uv_mod_cont (x::p)= gcd_int x (uv_mod_cont(p)); neuper@37950: neuper@37950: (*. divides each coefficient of a uv polynomial by y .*) neuper@37950: fun uv_mod_div_list (p:uv_poly,0) = raise error("UV_MOD_DIV_LIST_EXCEPTION: Division by zero") neuper@37950: | uv_mod_div_list ([],y) = []:uv_poly neuper@37950: | uv_mod_div_list (x::p,y) = (x div y)::uv_mod_div_list(p,y); neuper@37950: neuper@37950: (*. calculates the primitiv part of a uv polynomial .*) neuper@37950: fun uv_mod_pp ([]:uv_poly) = []:uv_poly neuper@37950: | uv_mod_pp p = neuper@37950: let neuper@38006: val c= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: c:=uv_mod_cont(p); neuper@37950: neuper@37950: if !c=0 then raise error ("RATIONALS_UV_MOD_PP_EXCEPTION: content is 0") neuper@37950: else uv_mod_div_list(p,!c) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. gets the leading coefficient of a uv polynomial .*) neuper@37950: fun uv_mod_lc ([]:uv_poly) = 0 neuper@37950: | uv_mod_lc p = hd(rev(p)); neuper@37950: neuper@37950: (*. calculates the euklidean polynomial remainder sequence in Zp .*) neuper@37950: fun uv_mod_prs_euklid_p(p1:uv_poly,p2:uv_poly,p)= neuper@37950: let neuper@38006: val f = Unsynchronized.ref []; neuper@38006: val f'= Unsynchronized.ref p2; neuper@38006: val fi= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: f:=p2::p1::[]; neuper@37950: while uv_mod_deg(!f')>0 do neuper@37950: ( neuper@37950: f':=uv_mod_prestp (hd(tl(!f))) (hd(!f)) p; neuper@37950: if (!f')<>[] then neuper@37950: ( neuper@37950: fi:=(!f'); neuper@37950: f:=(!fi)::(!f) neuper@37950: ) neuper@37950: else () neuper@37950: ); neuper@37950: (!f) neuper@37950: neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the gcd of p1 and p2 in Zp .*) neuper@37950: fun uv_mod_gcd_modp ([]:uv_poly) (p2:uv_poly) p = p2:uv_poly neuper@37950: | uv_mod_gcd_modp p1 [] p= p1 neuper@37950: | uv_mod_gcd_modp p1 p2 p= neuper@37950: let neuper@38006: val p1'= Unsynchronized.ref []; neuper@38006: val p2'= Unsynchronized.ref []; neuper@38006: val pc= Unsynchronized.ref []; neuper@38006: val g= Unsynchronized.ref []; neuper@38006: val d= Unsynchronized.ref 0; neuper@38006: val prs= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: if uv_mod_deg(p1)>=uv_mod_deg(p2) then neuper@37950: ( neuper@37950: p1':=uv_mod_list_modp (uv_mod_pp(p1)) p; neuper@37950: p2':=uv_mod_list_modp (uv_mod_pp(p2)) p neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=uv_mod_list_modp (uv_mod_pp(p2)) p; neuper@37950: p2':=uv_mod_list_modp (uv_mod_pp(p1)) p neuper@37950: ); neuper@37950: d:=uv_mod_mod2((gcd_int (uv_mod_cont(p1))) (uv_mod_cont(p2)), p) ; neuper@37950: if !d>(p div 2) then d:=(!d)-p else (); neuper@37950: neuper@37950: prs:=uv_mod_prs_euklid_p(!p1',!p2',p); neuper@37950: neuper@37950: if hd(!prs)=[] then pc:=hd(tl(!prs)) neuper@37950: else pc:=hd(!prs); neuper@37950: neuper@37950: g:=uv_mod_smul_poly(uv_mod_pp(!pc),!d); neuper@37950: !g neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the minimum of two real values x and y .*) neuper@37978: fun uv_mod_r_min(x,y):Real.real = if x>y then y else x; neuper@37950: neuper@37950: (*. calculates the minimum of two integer values x and y .*) neuper@37950: fun uv_mod_min(x,y) = if x>y then y else x; neuper@37950: neuper@37950: (*. adds the squared values of a integer list .*) neuper@37950: fun uv_mod_add_qu [] = 0.0 neuper@37978: | uv_mod_add_qu (x::p) = Real.fromInt(x)*Real.fromInt(x) + uv_mod_add_qu p; neuper@37950: neuper@37950: (*. calculates the euklidean norm .*) neuper@37950: fun uv_mod_norm ([]:uv_poly) = 0.0 neuper@37950: | uv_mod_norm p = Math.sqrt(uv_mod_add_qu(p)); neuper@37950: neuper@37950: (*. multipies two values a and b .*) neuper@37950: fun uv_mod_multi a b = a * b; neuper@37950: neuper@37950: (*. decides if x is a prim, the list contains all primes which are lower then x .*) neuper@37950: fun uv_mod_prim(x,[])= false neuper@37950: | uv_mod_prim(x,[y])=if ((x mod y) <> 0) then true neuper@37950: else false neuper@37950: | uv_mod_prim(x,y::ys) = if uv_mod_prim(x,[y]) neuper@37950: then neuper@37950: if uv_mod_prim(x,ys) then true neuper@37950: else false neuper@37950: else false; neuper@37950: neuper@37950: (*. gets the first prime, which is greater than p and does not divide g .*) neuper@37950: fun uv_mod_nextprime(g,p)= neuper@37950: let neuper@38006: val list= Unsynchronized.ref [2]; neuper@38006: val exit= Unsynchronized.ref 0; neuper@38006: val i = Unsynchronized.ref 2 neuper@37950: in neuper@37950: while (!i 0) neuper@37950: then neuper@37950: ( neuper@37950: list:= (!i)::(!list); neuper@37950: i:= (!i)+1 neuper@37950: ) neuper@37950: else i:=(!i)+1 neuper@37950: ) neuper@37950: else i:= (!i)+1 neuper@37950: ); neuper@37950: i:=(p+1); neuper@37950: while (!exit=0) do (* calculate next prime which does not divide g *) neuper@37950: ( neuper@37950: if uv_mod_prim(!i,!list) then neuper@37950: ( neuper@37950: if (g mod !i <> 0) neuper@37950: then neuper@37950: ( neuper@37950: list:= (!i)::(!list); neuper@37950: exit:= (!i) neuper@37950: ) neuper@37950: else i:=(!i)+1 neuper@37950: ) neuper@37950: else i:= (!i)+1 neuper@37950: ); neuper@37950: !exit neuper@37950: end; neuper@37950: neuper@37950: (*. decides if p1 is a factor of p2 in Zp .*) neuper@37950: fun uv_mod_dividesp ([]:uv_poly) (p2:uv_poly) p= raise error("UV_MOD_DIVIDESP: Division by zero") neuper@37950: | uv_mod_dividesp p1 p2 p= if uv_mod_prestp p2 p1 p = [] then true else false; neuper@37950: neuper@37950: (*. decides if p1 is a factor of p2 .*) neuper@37950: fun uv_mod_divides ([]:uv_poly) (p2:uv_poly) = raise error("UV_MOD_DIVIDES: Division by zero") neuper@37950: | uv_mod_divides p1 p2 = if uv_mod_prest p2 p1 = [] then true else false; neuper@37950: neuper@37950: (*. chinese remainder algorithm .*) neuper@37950: fun uv_mod_cra2(r1,r2,m1,m2)= neuper@37950: let neuper@38006: val c= Unsynchronized.ref 0; neuper@38006: val r1'= Unsynchronized.ref 0; neuper@38006: val d= Unsynchronized.ref 0; neuper@38006: val a= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: while (uv_mod_mod2((!c)*m1,m2))<>1 do neuper@37950: ( neuper@37950: c:=(!c)+1 neuper@37950: ); neuper@37950: r1':= uv_mod_mod2(r1,m1); neuper@37950: d:=uv_mod_mod2(((r2-(!r1'))*(!c)),m2); neuper@37950: !r1'+(!d)*m1 neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. applies the chinese remainder algorithmen to the coefficients of x1 and x2 .*) neuper@37950: fun uv_mod_cra_2 ([],[],m1,m2) = [] neuper@37950: | uv_mod_cra_2 ([],x2,m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x1") neuper@37950: | uv_mod_cra_2 (x1,[],m1,m2) = raise error("UV_MOD_CRA_2_EXCEPTION: invalid call x2") neuper@37950: | 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)); neuper@37950: neuper@37950: (*. calculates the gcd of two uv polynomials p1' and p2' with the modular algorithm .*) neuper@37950: fun uv_mod_gcd (p1':uv_poly) (p2':uv_poly) = neuper@37950: let neuper@38006: val p1= Unsynchronized.ref (uv_mod_pp(p1')); neuper@38006: val p2= Unsynchronized.ref (uv_mod_pp(p2')); neuper@37950: val c=gcd_int (uv_mod_cont(p1')) (uv_mod_cont(p2')); neuper@38006: val temp= Unsynchronized.ref []; neuper@38006: val cp= Unsynchronized.ref []; neuper@38006: val qp= Unsynchronized.ref []; neuper@38006: val q= Unsynchronized.ref []; neuper@38006: val pn= Unsynchronized.ref 0; neuper@38006: val d= Unsynchronized.ref 0; neuper@38006: val g1= Unsynchronized.ref 0; neuper@38006: val p= Unsynchronized.ref 0; neuper@38006: val m= Unsynchronized.ref 0; neuper@38006: val exit= Unsynchronized.ref 0; neuper@38006: val i= Unsynchronized.ref 1; neuper@37950: in neuper@37950: if length(!p1)>length(!p2) then () neuper@37950: else neuper@37950: ( neuper@37950: temp:= !p1; neuper@37950: p1:= !p2; neuper@37950: p2:= !temp neuper@37950: ); neuper@37950: neuper@37950: neuper@37950: d:=gcd_int (uv_mod_lc(!p1)) (uv_mod_lc(!p2)); neuper@37950: g1:=uv_mod_lc(!p1)*uv_mod_lc(!p2); neuper@37950: p:=4; neuper@37950: neuper@37978: m:=Real.ceil(2.0 * Real.fromInt(!d) * neuper@37978: Real.fromInt(power 2 (uv_mod_min(uv_mod_deg(!p2),uv_mod_deg(!p1)))) * neuper@37978: Real.fromInt(!d) * neuper@37978: uv_mod_r_min(uv_mod_norm(!p1) / Real.fromInt(abs(uv_mod_lc(!p1))), neuper@37978: uv_mod_norm(!p2) / Real.fromInt(abs(uv_mod_lc(!p2))))); neuper@37950: neuper@37950: while (!exit=0) do neuper@37950: ( neuper@37950: p:=uv_mod_nextprime(!d,!p); neuper@37950: cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)) ; neuper@37950: if abs(uv_mod_lc(!cp))<>1 then (* leading coefficient = 1 ? *) neuper@37950: ( neuper@37950: i:=1; neuper@37950: while (!i)<(!p) andalso (abs(uv_mod_mod2((uv_mod_lc(!cp)*(!i)),(!p)))<>1) do neuper@37950: ( neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p) neuper@37950: ) neuper@37950: else (); neuper@37950: neuper@37950: qp:= ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp)); neuper@37950: neuper@37950: if uv_mod_deg(!qp)=0 then (q:=[1]; exit:=1) else (); neuper@37950: neuper@37950: pn:=(!p); neuper@37950: q:=(!qp); neuper@37950: neuper@37950: while !pn<= !m andalso !m>(!p) andalso !exit=0 do neuper@37950: ( neuper@37950: p:=uv_mod_nextprime(!d,!p); neuper@37950: cp:=(uv_mod_gcd_modp (uv_mod_list_modp(!p1) (!p)) (uv_mod_list_modp(!p2) (!p)) (!p)); neuper@37950: if uv_mod_lc(!cp)<>1 then (* leading coefficient = 1 ? *) neuper@37950: ( neuper@37950: i:=1; neuper@37950: while (!i)<(!p) andalso ((uv_mod_mod2((uv_mod_lc(!q)*(!i)),(!p)))<>1) do neuper@37950: ( neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: cp:=uv_mod_list_modp (map (uv_mod_multi (!i)) (!cp)) (!p) neuper@37950: ) neuper@37950: else (); neuper@37950: neuper@37950: qp:=uv_mod_list_modp ((map (uv_mod_multi (uv_mod_mod2(!d,!p)))) (!cp) ) (!p); neuper@37950: if uv_mod_deg(!qp)>uv_mod_deg(!q) then neuper@37950: ( neuper@37950: (*print("degree to high!!!\n")*) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if uv_mod_deg(!qp)=uv_mod_deg(!q) then neuper@37950: ( neuper@37950: q:=uv_mod_cra_2(!q,!qp,!pn,!p); neuper@37950: pn:=(!pn) * !p; neuper@37950: q:=uv_mod_pp(uv_mod_list_modp (!q) (!pn)); (* found already gcd ? *) neuper@37950: if (uv_mod_divides (!q) (p1')) andalso (uv_mod_divides (!q) (p2')) then (exit:=1) else () neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if uv_mod_deg(!qp) [(5,[5,3,0]),(4,[3,0,2]),(2,[2,1,3]),(~1,[0,0,1]),(~19,[0,0,0])] .*) neuper@37950: neuper@37950: (*. global variables .*) neuper@37950: (*. order indicators .*) neuper@37950: val LEX_=0; (* lexicographical term order *) neuper@37950: val GGO_=1; (* greatest degree order *) neuper@37950: neuper@37950: (*. datatypes for internal representation.*) neuper@37950: type mv_monom = (int * (*.coefficient or the monom.*) neuper@37950: int list); (*.list of exponents) .*) neuper@37950: fun mv_monom2str (i, is) = "("^ int2str i^"," ^ ints2str' is ^ ")"; neuper@37950: neuper@37950: type mv_poly = mv_monom list; neuper@37950: fun mv_poly2str p = (strs2str' o (map mv_monom2str)) p; neuper@37950: neuper@37950: (*. help function for monom_greater and geq .*) neuper@37950: fun mv_mg_hlp([]) = EQUAL neuper@37950: | mv_mg_hlp(x::list)=if x<0 then LESS neuper@37950: else if x>0 then GREATER neuper@37950: else mv_mg_hlp(list); neuper@37950: neuper@37950: (*. adds a list of values .*) neuper@37950: fun mv_addlist([]) = 0 neuper@37950: | mv_addlist(p1) = hd(p1)+mv_addlist(tl(p1)); neuper@37950: neuper@37950: (*. tests if the monomial M1 is greater as the monomial M2 and returns a boolean value .*) neuper@37950: (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*) neuper@37950: fun mv_monom_greater((M1x,M1l):mv_monom,(M2x,M2l):mv_monom,order)= neuper@37950: if order=LEX_ then neuper@37950: ( neuper@37950: if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error") neuper@37950: else if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true neuper@37950: ) neuper@37950: else neuper@37950: if order=GGO_ then neuper@37950: ( neuper@37950: if length(M1l)<>length(M2l) then raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Order error") neuper@37950: else neuper@37950: if mv_addlist(M1l)=mv_addlist(M2l) then if (mv_mg_hlp((map op- (M1l~~M2l)))<>GREATER) then false else true neuper@37950: else if mv_addlist(M1l)>mv_addlist(M2l) then true else false neuper@37950: ) neuper@37950: else raise error ("RATIONALS_MV_MONOM_GREATER_EXCEPTION: Wrong Order"); neuper@37950: neuper@37950: (*. tests if the monomial X is greater as the monomial Y and returns a order value (GREATER,EQUAL,LESS) .*) neuper@37950: (*. 2 orders are implemented LEX_/GGO_ (lexigraphical/greatest degree order) .*) neuper@37950: fun mv_geq order ((x1,x):mv_monom,(x2,y):mv_monom) = neuper@37950: let neuper@38006: val temp= Unsynchronized.ref EQUAL; neuper@37950: in neuper@37950: if order=LEX_ then neuper@37950: ( neuper@37950: if length(x)<>length(y) then neuper@37950: raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error") neuper@37950: else neuper@37950: ( neuper@37950: temp:=mv_mg_hlp((map op- (x~~y))); neuper@37950: if !temp=EQUAL then neuper@37950: ( if x1=x2 then EQUAL neuper@37950: else if x1>x2 then GREATER neuper@37950: else LESS neuper@37950: ) neuper@37950: else (!temp) neuper@37950: ) neuper@37950: ) neuper@37950: else neuper@37950: if order=GGO_ then neuper@37950: ( neuper@37950: if length(x)<>length(y) then neuper@37950: raise error ("RATIONALS_MV_GEQ_EXCEPTION: Order error") neuper@37950: else neuper@37950: if mv_addlist(x)=mv_addlist(y) then neuper@37950: (mv_mg_hlp((map op- (x~~y)))) neuper@37950: else if mv_addlist(x)>mv_addlist(y) then GREATER else LESS neuper@37950: ) neuper@37950: else raise error ("RATIONALS_MV_GEQ_EXCEPTION: Wrong Order") neuper@37950: end; neuper@37950: neuper@37950: (*. cuts the first variable from a polynomial .*) neuper@37950: fun mv_cut([]:mv_poly)=[]:mv_poly neuper@37950: | mv_cut((x,[])::list) = raise error ("RATIONALS_MV_CUT_EXCEPTION: Invalid list ") neuper@37950: | mv_cut((x,y::ys)::list)=(x,ys)::mv_cut(list); neuper@37950: neuper@37950: (*. leading power product .*) neuper@37950: fun mv_lpp([]:mv_poly,order) = [] neuper@37950: | mv_lpp([(x,y)],order) = y neuper@37950: | mv_lpp(p1,order) = #2(hd(rev(sort (mv_geq order) p1))); neuper@37950: neuper@37950: (*. leading monomial .*) neuper@37950: fun mv_lm([]:mv_poly,order) = (0,[]):mv_monom neuper@37950: | mv_lm([x],order) = x neuper@37950: | mv_lm(p1,order) = hd(rev(sort (mv_geq order) p1)); neuper@37950: neuper@37950: (*. leading coefficient in term order .*) neuper@37950: fun mv_lc2([]:mv_poly,order) = 0 neuper@37950: | mv_lc2([(x,y)],order) = x neuper@37950: | mv_lc2(p1,order) = #1(hd(rev(sort (mv_geq order) p1))); neuper@37950: neuper@37950: neuper@37950: (*. reverse the coefficients in mv polynomial .*) neuper@37950: fun mv_rev_to([]:mv_poly) = []:mv_poly neuper@37950: | mv_rev_to((c,e)::xs) = (c,rev(e))::mv_rev_to(xs); neuper@37950: neuper@37950: (*. leading coefficient in reverse term order .*) neuper@37950: fun mv_lc([]:mv_poly,order) = []:mv_poly neuper@37950: | mv_lc([(x,y)],order) = mv_rev_to(mv_cut(mv_rev_to([(x,y)]))) neuper@37950: | mv_lc(p1,order) = neuper@37950: let neuper@38006: val p1o= Unsynchronized.ref (rev(sort (mv_geq order) (mv_rev_to(p1)))); neuper@37950: val lp=hd(#2(hd(!p1o))); neuper@38006: val lc= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: while (length(!p1o)>0 andalso hd(#2(hd(!p1o)))=lp) do neuper@37950: ( neuper@37950: lc:=hd(mv_cut([hd(!p1o)]))::(!lc); neuper@37950: p1o:=tl(!p1o) neuper@37950: ); neuper@37950: if !lc=[] then raise error ("RATIONALS_MV_LC_EXCEPTION: lc is empty") else (); neuper@37950: mv_rev_to(!lc) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. compares two powerproducts .*) neuper@37950: fun mv_monom_equal((_,xlist):mv_monom,(_,ylist):mv_monom) = (foldr and_) (((map op=) (xlist~~ylist)),true); neuper@37950: neuper@37950: (*. help function for mv_add .*) neuper@37950: fun mv_madd([]:mv_poly,[]:mv_poly,order) = []:mv_poly neuper@37950: | mv_madd([(0,_)],p2,order) = p2 neuper@37950: | mv_madd(p1,[(0,_)],order) = p1 neuper@37950: | mv_madd([],p2,order) = p2 neuper@37950: | mv_madd(p1,[],order) = p1 neuper@37950: | mv_madd(p1,p2,order) = neuper@37950: ( neuper@37950: if mv_monom_greater(hd(p1),hd(p2),order) neuper@37950: then hd(p1)::mv_madd(tl(p1),p2,order) neuper@37950: else if mv_monom_equal(hd(p1),hd(p2)) neuper@37950: then if mv_lc2(p1,order)+mv_lc2(p2,order)<>0 neuper@37950: then (mv_lc2(p1,order)+mv_lc2(p2,order),mv_lpp(p1,order))::mv_madd(tl(p1),tl(p2),order) neuper@37950: else mv_madd(tl(p1),tl(p2),order) neuper@37950: else hd(p2)::mv_madd(p1,tl(p2),order) neuper@37950: ) neuper@37950: neuper@37950: (*. adds two multivariate polynomials .*) neuper@37950: fun mv_add([]:mv_poly,p2:mv_poly,order) = p2 neuper@37950: | mv_add(p1,[],order) = p1 neuper@37950: | mv_add(p1,p2,order) = mv_madd(rev(sort (mv_geq order) p1),rev(sort (mv_geq order) p2), order); neuper@37950: neuper@37950: (*. monom multiplication .*) neuper@37950: fun mv_mmul((x1,y1):mv_monom,(x2,y2):mv_monom)=(x1*x2,(map op+) (y1~~y2)):mv_monom; neuper@37950: neuper@37950: (*. deletes all monomials with coefficient 0 .*) neuper@37950: fun mv_shorten([]:mv_poly,order) = []:mv_poly neuper@37950: | mv_shorten(x::xs,order)=mv_madd([x],mv_shorten(xs,order),order); neuper@37950: neuper@37950: (*. zeros a list .*) neuper@37950: fun mv_null2([])=[] neuper@37950: | mv_null2(x::l)=0::mv_null2(l); neuper@37950: neuper@37950: (*. multiplies two multivariate polynomials .*) neuper@37950: fun mv_mul([]:mv_poly,[]:mv_poly,_) = []:mv_poly neuper@37950: | mv_mul([],y::p2,_) = [(0,mv_null2(#2(y)))] neuper@37950: | mv_mul(x::p1,[],_) = [(0,mv_null2(#2(x)))] neuper@37950: | mv_mul(x::p1,y::p2,order) = mv_shorten(rev(sort (mv_geq order) (mv_mmul(x,y) :: (mv_mul(p1,y::p2,order) @ neuper@37950: mv_mul([x],p2,order)))),order); neuper@37950: neuper@37950: (*. gets the maximum value of a list .*) neuper@37950: fun mv_getmax([])=0 neuper@37950: | mv_getmax(x::p1)= let neuper@37950: val m=mv_getmax(p1); neuper@37950: in neuper@37950: if m>x then m neuper@37950: else x neuper@37950: end; neuper@37950: (*. calculates the maximum degree of an multivariate polynomial .*) neuper@37950: fun mv_grad([]:mv_poly) = 0 neuper@37950: | mv_grad(p1:mv_poly)= mv_getmax((map mv_addlist) ((map #2) p1)); neuper@37950: neuper@37950: (*. converts the sign of a value .*) neuper@37950: fun mv_minus(x)=(~1) * x; neuper@37950: neuper@37950: (*. converts the sign of all coefficients of a polynomial .*) neuper@37950: fun mv_minus2([]:mv_poly)=[]:mv_poly neuper@37950: | mv_minus2(p1)=(mv_minus(#1(hd(p1))),#2(hd(p1)))::(mv_minus2(tl(p1))); neuper@37950: neuper@37950: (*. searches for a negativ value in a list .*) neuper@37950: fun mv_is_negativ([])=false neuper@37950: | mv_is_negativ(x::xs)=if x<0 then true else mv_is_negativ(xs); neuper@37950: neuper@37950: (*. division of monomials .*) neuper@37950: fun mv_mdiv((0,[]):mv_monom,_:mv_monom)=(0,[]):mv_monom neuper@37950: | mv_mdiv(_,(0,[]))= raise error ("RATIONALS_MV_MDIV_EXCEPTION Division by 0 ") neuper@37950: | mv_mdiv(p1:mv_monom,p2:mv_monom)= neuper@37950: let neuper@38006: val c= Unsynchronized.ref (#1(p2)); neuper@38006: val pp= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: if !c=0 then raise error("MV_MDIV_EXCEPTION Dividing by zero") neuper@37950: else c:=(#1(p1) div #1(p2)); neuper@37950: if #1(p2)<>0 then neuper@37950: ( neuper@37950: pp:=(#2(mv_mmul((1,#2(p1)),(1,(map mv_minus) (#2(p2)))))); neuper@37950: if mv_is_negativ(!pp) then (0,!pp) neuper@37950: else (!c,!pp) neuper@37950: ) neuper@37950: else raise error("MV_MDIV_EXCEPTION Dividing by empty Polynom") neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. prints a polynom for (internal use only) .*) neuper@37978: fun mv_print_poly([]:mv_poly)=writeln("[]\n") neuper@37978: | mv_print_poly((x,y)::[])= writeln("("^Int.toString(x)^","^ints2str(y)^")\n") neuper@37978: | mv_print_poly((x,y)::p1) = (writeln("("^Int.toString(x)^","^ints2str(y)^"),");mv_print_poly(p1)); neuper@37950: neuper@37950: neuper@37950: (*. division of two multivariate polynomials .*) neuper@37950: fun mv_division([]:mv_poly,g:mv_poly,order)=([]:mv_poly,[]:mv_poly) neuper@37950: | mv_division(f,[],order)= raise error ("RATIONALS_MV_DIVISION_EXCEPTION Division by zero") neuper@37950: | mv_division(f,g,order)= neuper@37950: let neuper@38006: val r= Unsynchronized.ref []; neuper@38006: val q= Unsynchronized.ref []; neuper@38006: val g'= Unsynchronized.ref ([] : mv_monom list); neuper@38006: val k= Unsynchronized.ref 0; neuper@38006: val m= Unsynchronized.ref (0,[0]); neuper@38006: val exit= Unsynchronized.ref 0; neuper@37950: in neuper@37950: r := rev(sort (mv_geq order) (mv_shorten(f,order))); neuper@37950: g':= rev(sort (mv_geq order) (mv_shorten(g,order))); neuper@37950: if #1(hd(!g'))=0 then raise error("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero") else (); neuper@37950: if (mv_monom_greater (hd(!g'),hd(!r),order)) then ([(0,mv_null2(#2(hd(f))))],(!r)) neuper@37950: else neuper@37950: ( neuper@37950: exit:=0; neuper@37950: while (if (!exit)=0 then not(mv_monom_greater (hd(!g'),hd(!r),order)) else false) do neuper@37950: ( neuper@37950: if (#1(mv_lm(!g',order)))<>0 then m:=mv_mdiv(mv_lm(!r,order),mv_lm(!g',order)) neuper@37950: else raise error ("RATIONALS_MV_DIVISION_EXCEPTION: Dividing by zero"); neuper@37950: if #1(!m)<>0 then neuper@37950: ( neuper@37950: q:=(!m)::(!q); neuper@37950: r:=mv_add((!r),mv_minus2(mv_mul(!g',[!m],order)),order) neuper@37950: ) neuper@37950: else exit:=1; neuper@37950: if (if length(!r)<>0 then length(!g')<>0 else false) then () neuper@37950: else (exit:=1) neuper@37950: ); neuper@37950: (rev(!q),!r) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. multiplies a polynomial with an integer .*) neuper@37950: fun mv_skalar_mul([]:mv_poly,c) = []:mv_poly neuper@37950: | mv_skalar_mul((x,y)::p1,c) = ((x * c),y)::mv_skalar_mul(p1,c); neuper@37950: neuper@37950: (*. inserts the a first variable into an polynomial with exponent v .*) neuper@37950: fun mv_correct([]:mv_poly,v:int)=[]:mv_poly neuper@37950: | mv_correct((x,y)::list,v:int)=(x,v::y)::mv_correct(list,v); neuper@37950: neuper@37950: (*. multivariate case .*) neuper@37950: neuper@37950: (*. decides if x is a factor of y .*) neuper@37950: fun mv_divides([]:mv_poly,[]:mv_poly)= raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero") neuper@37950: | mv_divides(x,[]) = raise error("RATIONALS_MV_DIVIDES_EXCEPTION: division by zero") neuper@37950: | mv_divides(x:mv_poly,y:mv_poly) = #2(mv_division(y,x,LEX_))=[]; neuper@37950: neuper@37950: (*. gets the maximum of a and b .*) neuper@37950: fun mv_max(a,b) = if a>b then a else b; neuper@37950: neuper@37950: (*. gets the maximum exponent of a mv polynomial in the lexicographic term order .*) neuper@37950: fun mv_deg([]:mv_poly) = 0 neuper@37950: | mv_deg(p1)= neuper@37950: let neuper@37950: val p1'=mv_shorten(p1,LEX_); neuper@37950: in neuper@37950: if length(p1')=0 then 0 neuper@37950: else mv_max(hd(#2(hd(p1'))),mv_deg(tl(p1'))) neuper@37950: end; neuper@37950: neuper@37950: (*. gets the maximum exponent of a mv polynomial in the reverse lexicographic term order .*) neuper@37950: fun mv_deg2([]:mv_poly) = 0 neuper@37950: | mv_deg2(p1)= neuper@37950: let neuper@37950: val p1'=mv_shorten(p1,LEX_); neuper@37950: in neuper@37950: if length(p1')=0 then 0 neuper@37950: else mv_max(hd(rev(#2(hd(p1')))),mv_deg2(tl(p1'))) neuper@37950: end; neuper@37950: neuper@37950: (*. evaluates the mv polynomial at the value v of the main variable .*) neuper@37950: fun mv_subs([]:mv_poly,v) = []:mv_poly neuper@37950: | mv_subs((c,e)::p1:mv_poly,v) = mv_skalar_mul(mv_cut([(c,e)]),power v (hd(e))) @ mv_subs(p1,v); neuper@37950: neuper@37950: (*. calculates the content of a uv-polynomial in mv-representation .*) neuper@37950: fun uv_content2([]:mv_poly) = 0 neuper@37950: | uv_content2((c,e)::p1) = (gcd_int c (uv_content2(p1))); neuper@37950: neuper@37950: (*. converts a uv-polynomial from mv-representation to uv-representation .*) neuper@37950: fun uv_to_list ([]:mv_poly)=[]:uv_poly neuper@37950: | uv_to_list ((c1,e1)::others) = neuper@37950: let neuper@38006: val count= Unsynchronized.ref 0; neuper@37950: val max=mv_grad((c1,e1)::others); neuper@38006: val help= Unsynchronized.ref ((c1,e1)::others); neuper@38006: val list= Unsynchronized.ref []; neuper@37950: in neuper@37950: if length(e1)>1 then raise error ("RATIONALS_TO_LIST_EXCEPTION: not univariate") neuper@37950: else if length(e1)=0 then [c1] neuper@37950: else neuper@37950: ( neuper@37950: count:=0; neuper@37950: while (!count)<=max do neuper@37950: ( neuper@37950: if length(!help)>0 andalso hd(#2(hd(!help)))=max-(!count) then neuper@37950: ( neuper@37950: list:=(#1(hd(!help)))::(!list); neuper@37950: help:=tl(!help) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: list:= 0::(!list) neuper@37950: ); neuper@37950: count := (!count) + 1 neuper@37950: ); neuper@37950: (!list) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. converts a uv-polynomial from uv-representation to mv-representation .*) neuper@37950: fun uv_to_poly ([]:uv_poly) = []:mv_poly neuper@37950: | uv_to_poly p1 = neuper@37950: let neuper@38006: val count= Unsynchronized.ref 0; neuper@38006: val help= Unsynchronized.ref p1; neuper@38006: val list= Unsynchronized.ref []; neuper@37950: in neuper@37950: while length(!help)>0 do neuper@37950: ( neuper@37950: if hd(!help)=0 then () neuper@37950: else list:=(hd(!help),[!count])::(!list); neuper@37950: count:=(!count)+1; neuper@37950: help:=tl(!help) neuper@37950: ); neuper@37950: (!list) neuper@37950: end; neuper@37950: neuper@37950: (*. univariate gcd calculation from polynomials in multivariate representation .*) neuper@37950: fun uv_gcd ([]:mv_poly) p2 = p2 neuper@37950: | uv_gcd p1 ([]:mv_poly) = p1 neuper@37950: | uv_gcd p1 [(c,[e])] = neuper@37950: let neuper@38006: val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_)))); neuper@37950: val min=uv_mod_min(e,(hd(#2(hd(rev(!list)))))); neuper@37950: in neuper@37950: [(gcd_int (uv_content2(p1)) c,[min])] neuper@37950: end neuper@37950: | uv_gcd [(c,[e])] p2 = neuper@37950: let neuper@38006: val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p2,LEX_)))); neuper@37950: val min=uv_mod_min(e,(hd(#2(hd(rev(!list)))))); neuper@37950: in neuper@37950: [(gcd_int (uv_content2(p2)) c,[min])] neuper@37950: end neuper@37950: | uv_gcd p11 p22 = uv_to_poly(uv_mod_gcd (uv_to_list(mv_shorten(p11,LEX_))) (uv_to_list(mv_shorten(p22,LEX_)))); neuper@37950: neuper@37950: (*. help function for the newton interpolation .*) neuper@37950: fun mv_newton_help ([]:mv_poly list,k:int) = []:mv_poly list neuper@37950: | mv_newton_help (pl:mv_poly list,k) = neuper@37950: let neuper@38006: val x= Unsynchronized.ref (rev(pl)); neuper@38006: val t= Unsynchronized.ref []; neuper@38006: val y= Unsynchronized.ref []; neuper@38006: val n= Unsynchronized.ref 1; neuper@38006: val n1= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: while length(!x)>1 do neuper@37950: ( neuper@37950: if length(hd(!x))>0 then n1:=mv_null2(#2(hd(hd(!x)))) neuper@37950: else if length(hd(tl(!x)))>0 then n1:=mv_null2(#2(hd(hd(tl(!x))))) neuper@37950: else n1:=[]; neuper@37950: t:= #1(mv_division(mv_add(hd(!x),mv_skalar_mul(hd(tl(!x)),~1),LEX_),[(k,!n1)],LEX_)); neuper@37950: y:=(!t)::(!y); neuper@37950: x:=tl(!x) neuper@37950: ); neuper@37950: (!y) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. help function for the newton interpolation .*) neuper@37950: fun mv_newton_add ([]:mv_poly list) t = []:mv_poly neuper@37950: | mv_newton_add [x:mv_poly] t = x neuper@37950: | mv_newton_add (pl:mv_poly list) t = neuper@37950: let neuper@38006: val expos= Unsynchronized.ref []; neuper@38006: val pll= Unsynchronized.ref pl; neuper@37950: in neuper@37950: ( neuper@37950: neuper@37950: while length(!pll)>0 andalso hd(!pll)=[] do neuper@37950: ( neuper@37950: pll:=tl(!pll) neuper@37950: ); neuper@37950: if length(!pll)>0 then expos:= #2(hd(hd(!pll))) else expos:=[]; neuper@37950: mv_add(hd(pl), neuper@37950: mv_mul( neuper@37950: mv_add(mv_correct(mv_cut([(1,mv_null2(!expos))]),1),[(~t,mv_null2(!expos))],LEX_), neuper@37950: mv_newton_add (tl(pl)) (t+1), neuper@37950: LEX_ neuper@37950: ), neuper@37950: LEX_) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the newton interpolation with polynomial coefficients .*) neuper@37950: (*. step-depth is 1 and if the result is not an integerpolynomial .*) neuper@37950: (*. this function returns [] .*) neuper@37950: fun mv_newton ([]:(mv_poly) list) = []:mv_poly neuper@37950: | mv_newton ([mp]:(mv_poly) list) = mp:mv_poly neuper@37950: | mv_newton pl = neuper@37950: let neuper@38006: val c= Unsynchronized.ref pl; neuper@38006: val c1= Unsynchronized.ref []; neuper@37950: val n=length(pl); neuper@38006: val k= Unsynchronized.ref 1; neuper@38006: val i= Unsynchronized.ref n; neuper@38006: val ppl= Unsynchronized.ref []; neuper@37950: in neuper@37950: c1:=hd(pl)::[]; neuper@37950: c:=mv_newton_help(!c,!k); neuper@37950: c1:=(hd(!c))::(!c1); neuper@37950: while(length(!c)>1 andalso !k0 andalso hd(!c)=[] do c:=tl(!c); neuper@37950: if !c=[] then () else c:=mv_newton_help(!c,!k); neuper@37950: ppl:= !c; neuper@37950: if !c=[] then () else c1:=(hd(!c))::(!c1) neuper@37950: ); neuper@37950: while hd(!c1)=[] do c1:=tl(!c1); neuper@37950: c1:=rev(!c1); neuper@37950: ppl:= !c1; neuper@37950: mv_newton_add (!c1) 1 neuper@37950: end; neuper@37950: neuper@37950: (*. sets the exponents of the first variable to zero .*) neuper@37950: fun mv_null3([]:mv_poly) = []:mv_poly neuper@37950: | mv_null3((x,y)::xs) = (x,0::tl(y))::mv_null3(xs); neuper@37950: neuper@37950: (*. calculates the minimum exponents of a multivariate polynomial .*) neuper@37950: fun mv_min_pp([]:mv_poly)=[] neuper@37950: | mv_min_pp((c,e)::xs)= neuper@37950: let neuper@38006: val y= Unsynchronized.ref xs; neuper@38006: val x= Unsynchronized.ref []; neuper@37950: in neuper@37950: ( neuper@37950: x:=e; neuper@37950: while length(!y)>0 do neuper@37950: ( neuper@37950: x:=(map uv_mod_min) ((!x) ~~ (#2(hd(!y)))); neuper@37950: y:=tl(!y) neuper@37950: ); neuper@37950: !x neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. checks if all elements of the list have value zero .*) neuper@37950: fun list_is_null [] = true neuper@37950: | list_is_null (x::xs) = (x=0 andalso list_is_null(xs)); neuper@37950: neuper@37950: (* check if main variable is zero*) neuper@37950: fun main_zero (ms : mv_poly) = (list_is_null o (map (hd o #2))) ms; neuper@37950: neuper@37950: (*. calculates the content of an polynomial .*) neuper@37950: fun mv_content([]:mv_poly) = []:mv_poly neuper@37950: | mv_content(p1) = neuper@37950: let neuper@38006: val list= Unsynchronized.ref (rev(sort (mv_geq LEX_) (mv_shorten(p1,LEX_)))); neuper@38006: val test= Unsynchronized.ref (hd(#2(hd(!list)))); neuper@38006: val result= Unsynchronized.ref []; neuper@37950: val min=(hd(#2(hd(rev(!list))))); neuper@37950: in neuper@37950: ( neuper@37950: if length(!list)>1 then neuper@37950: ( neuper@37950: while (if length(!list)>0 then (hd(#2(hd(!list)))=(!test)) else false) do neuper@37950: ( neuper@37950: result:=(#1(hd(!list)),tl(#2(hd(!list))))::(!result); neuper@37950: neuper@37950: if length(!list)<1 then list:=[] neuper@37950: else list:=tl(!list) neuper@37950: neuper@37950: ); neuper@37950: if length(!list)>0 then neuper@37950: ( neuper@37950: list:=mv_gcd (!result) (mv_cut(mv_content(!list))) neuper@37950: ) neuper@37950: else list:=(!result); neuper@37950: list:=mv_correct(!list,0); neuper@37950: (!list) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: mv_null3(!list) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: neuper@37950: (*. calculates the primitiv part of a polynomial .*) neuper@37950: and mv_pp([]:mv_poly) = []:mv_poly neuper@37950: | mv_pp(p1) = let neuper@38006: val cont= Unsynchronized.ref []; neuper@38006: val pp= Unsynchronized.ref []; neuper@37950: in neuper@37950: cont:=mv_content(p1); neuper@37950: pp:=(#1(mv_division(p1,!cont,LEX_))); neuper@37950: if !pp=[] neuper@37950: then raise error("RATIONALS_MV_PP_EXCEPTION: Invalid Content ") neuper@37950: else (!pp) neuper@37950: end neuper@37950: neuper@37950: (*. calculates the gcd of two multivariate polynomials with a modular approach .*) neuper@37950: and mv_gcd ([]:mv_poly) ([]:mv_poly) :mv_poly= []:mv_poly neuper@37950: | mv_gcd ([]:mv_poly) (p2) :mv_poly= p2:mv_poly neuper@37950: | mv_gcd (p1:mv_poly) ([]) :mv_poly= p1:mv_poly neuper@37950: | mv_gcd ([(x,xs)]:mv_poly) ([(y,ys)]):mv_poly = neuper@37950: let neuper@37950: val xpoly:mv_poly = [(x,xs)]; neuper@37950: val ypoly:mv_poly = [(y,ys)]; neuper@37950: in neuper@37950: ( neuper@37950: if xs=ys then [((gcd_int x y),xs)] neuper@37950: else [((gcd_int x y),(map uv_mod_min)(xs~~ys))]:mv_poly neuper@37950: ) neuper@37950: end neuper@37950: | mv_gcd (p1:mv_poly) ([(y,ys)]) :mv_poly= neuper@37950: ( neuper@37950: [(gcd_int (uv_content2(p1)) (y),(map uv_mod_min)(mv_min_pp(p1)~~ys))]:mv_poly neuper@37950: ) neuper@37950: | mv_gcd ([(y,ys)]:mv_poly) (p2):mv_poly = neuper@37950: ( neuper@37950: [(gcd_int (uv_content2(p2)) (y),(map uv_mod_min)(mv_min_pp(p2)~~ys))]:mv_poly neuper@37950: ) neuper@37950: | mv_gcd (p1':mv_poly) (p2':mv_poly):mv_poly= neuper@37950: let neuper@37950: val vc=length(#2(hd(p1'))); neuper@37950: val cont = neuper@37950: ( neuper@37950: if main_zero(mv_content(p1')) andalso neuper@37950: (main_zero(mv_content(p2'))) then neuper@37950: mv_correct((mv_gcd (mv_cut(mv_content(p1'))) (mv_cut(mv_content(p2')))),0) neuper@37950: else neuper@37950: mv_gcd (mv_content(p1')) (mv_content(p2')) neuper@37950: ); neuper@37950: val p1= #1(mv_division(p1',mv_content(p1'),LEX_)); neuper@37950: val p2= #1(mv_division(p2',mv_content(p2'),LEX_)); neuper@38006: val gcd= Unsynchronized.ref []; neuper@38006: val candidate= Unsynchronized.ref []; neuper@38006: val interpolation_list= Unsynchronized.ref []; neuper@38006: val delta= Unsynchronized.ref []; neuper@38006: val p1r = Unsynchronized.ref []; neuper@38006: val p2r = Unsynchronized.ref []; neuper@38006: val p1r' = Unsynchronized.ref []; neuper@38006: val p2r' = Unsynchronized.ref []; neuper@38006: val factor= Unsynchronized.ref []; neuper@38006: val r= Unsynchronized.ref 0; neuper@38006: val gcd_r= Unsynchronized.ref []; neuper@38006: val d= Unsynchronized.ref 0; neuper@38006: val exit= Unsynchronized.ref 0; neuper@38006: val current_degree= Unsynchronized.ref 99999; (*. FIXME: unlimited ! .*) neuper@37950: in neuper@37950: ( neuper@37950: if vc<2 then (* areUnivariate(p1',p2') *) neuper@37950: ( neuper@37950: gcd:=uv_gcd (mv_shorten(p1',LEX_)) (mv_shorten(p2',LEX_)) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: while !exit=0 do neuper@37950: ( neuper@37950: r:=(!r)+1; neuper@37950: p1r := mv_lc(p1,LEX_); neuper@37950: p2r := mv_lc(p2,LEX_); neuper@37950: if main_zero(!p1r) andalso neuper@37950: main_zero(!p2r) neuper@37950: then neuper@37950: ( neuper@37950: delta := mv_correct((mv_gcd (mv_cut (!p1r)) (mv_cut (!p2r))),0) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: delta := mv_gcd (!p1r) (!p2r) neuper@37950: ); neuper@37950: (*if mv_shorten(mv_subs(!p1r,!r),LEX_)=[] andalso neuper@37950: mv_shorten(mv_subs(!p2r,!r),LEX_)=[] *) neuper@37950: if mv_lc2(mv_shorten(mv_subs(!p1r,!r),LEX_),LEX_)=0 andalso neuper@37950: mv_lc2(mv_shorten(mv_subs(!p2r,!r),LEX_),LEX_)=0 neuper@37950: then neuper@37950: ( neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: gcd_r:=mv_shorten(mv_gcd (mv_shorten(mv_subs(p1,!r),LEX_)) neuper@37950: (mv_shorten(mv_subs(p2,!r),LEX_)) ,LEX_); neuper@37950: gcd_r:= #1(mv_division(mv_mul(mv_correct(mv_subs(!delta,!r),0),!gcd_r,LEX_), neuper@37950: mv_correct(mv_lc(!gcd_r,LEX_),0),LEX_)); neuper@37950: d:=mv_deg2(!gcd_r); (* deg(gcd_r,z) *) neuper@37950: if (!d < !current_degree) then neuper@37950: ( neuper@37950: current_degree:= !d; neuper@37950: interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if (!d = !current_degree) then neuper@37950: ( neuper@37950: interpolation_list:=mv_correct(!gcd_r,0)::(!interpolation_list) neuper@37950: ) neuper@37950: else () neuper@37950: ) neuper@37950: ); neuper@37950: if length(!interpolation_list)> uv_mod_min(mv_deg(p1),mv_deg(p2)) then neuper@37950: ( neuper@37950: candidate := mv_newton(rev(!interpolation_list)); neuper@37950: if !candidate=[] then () neuper@37950: else neuper@37950: ( neuper@37950: candidate:=mv_pp(!candidate); neuper@37950: if mv_divides(!candidate,p1) andalso mv_divides(!candidate,p2) then neuper@37950: ( neuper@37950: gcd:= mv_mul(!candidate,cont,LEX_); neuper@37950: exit:=1 neuper@37950: ) neuper@37950: else () neuper@37950: ); neuper@37950: interpolation_list:=[mv_correct(!gcd_r,0)] neuper@37950: ) neuper@37950: else () neuper@37950: ) neuper@37950: ); neuper@37950: (!gcd):mv_poly neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: neuper@37950: (*. calculates the least common divisor of two polynomials .*) neuper@37950: fun mv_lcm (p1:mv_poly) (p2:mv_poly) :mv_poly = neuper@37950: ( neuper@37950: #1(mv_division(mv_mul(p1,p2,LEX_),mv_gcd p1 p2,LEX_)) neuper@37950: ); neuper@37950: neuper@37950: (*. gets the variables (strings) of a term .*) neuper@37950: fun get_vars(term1) = (map free2str) (vars term1); (*["a","b","c"]; *) neuper@37950: neuper@37950: (*. counts the negative coefficents in a polynomial .*) neuper@37950: fun count_neg ([]:mv_poly) = 0 neuper@37950: | count_neg ((c,e)::xs) = if c<0 then 1+count_neg xs neuper@37950: else count_neg xs; neuper@37950: neuper@37950: (*. help function for is_polynomial neuper@37950: checks the order of the operators .*) neuper@37950: fun test_polynomial (Const ("uminus",_) $ Free (str,_)) _ = true (*WN.13.3.03*) neuper@37950: | test_polynomial (t as Free(str,_)) v = true neuper@37950: | test_polynomial (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false neuper@37950: else (test_polynomial t1 "*") andalso (test_polynomial t2 "*") neuper@38014: | test_polynomial (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false neuper@37950: else (test_polynomial t1 " ") andalso (test_polynomial t2 " ") neuper@37950: | test_polynomial (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_polynomial t1 "^") andalso (test_polynomial t2 "^") neuper@37950: | test_polynomial _ v = false; neuper@37950: neuper@37950: (*. tests if a term is a polynomial .*) neuper@37950: fun is_polynomial t = test_polynomial t " "; neuper@37950: neuper@37950: (*. help function for is_expanded neuper@37950: checks the order of the operators .*) neuper@37950: fun test_exp (t as Free(str,_)) v = true neuper@37950: | test_exp (t as Const ("op *",_) $ t1 $ t2) v = if v="^" then false neuper@37950: else (test_exp t1 "*") andalso (test_exp t2 "*") neuper@38014: | test_exp (t as Const ("Groups.plus_class.plus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false neuper@37950: else (test_exp t1 " ") andalso (test_exp t2 " ") neuper@38014: | test_exp (t as Const ("Groups.minus_class.minus",_) $ t1 $ t2) v = if v="*" orelse v="^" then false neuper@37950: else (test_exp t1 " ") andalso (test_exp t2 " ") neuper@37950: | test_exp (t as Const ("Atools.pow",_) $ t1 $ t2) v = (test_exp t1 "^") andalso (test_exp t2 "^") neuper@37950: | test_exp _ v = false; neuper@37950: neuper@37950: neuper@37950: (*. help function for check_coeff: neuper@37950: converts the term to a list of coefficients .*) neuper@37950: fun term2coef' (t as Free(str,_(*typ*))) v :mv_poly option = neuper@37950: let neuper@38006: val x= Unsynchronized.ref NONE; neuper@38006: val len= Unsynchronized.ref 0; neuper@38006: val vl= Unsynchronized.ref []; neuper@38006: val vh= Unsynchronized.ref []; neuper@38006: val i= Unsynchronized.ref 0; neuper@37950: in neuper@37950: if is_numeral str then neuper@37950: ( neuper@37950: SOME [(((the o int_of_str) str),mv_null2(v))] handle _ => NONE neuper@37950: ) neuper@37950: else (* variable *) neuper@37950: ( neuper@37950: len:=length(v); neuper@37950: vh:=v; neuper@37950: while ((!len)>(!i)) do neuper@37950: ( neuper@37950: if str=hd((!vh)) then neuper@37950: ( neuper@37950: vl:=1::(!vl) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: vl:=0::(!vl) neuper@37950: ); neuper@37950: vh:=tl(!vh); neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: SOME [(1,rev(!vl))] handle _ => NONE neuper@37950: ) neuper@37950: end neuper@37950: | term2coef' (Const ("op *",_) $ t1 $ t2) v :mv_poly option= neuper@37950: let neuper@38006: val t1pp= Unsynchronized.ref []; neuper@38006: val t2pp= Unsynchronized.ref []; neuper@38006: val t1c= Unsynchronized.ref 0; neuper@38006: val t2c= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: t1pp:=(#2(hd(the(term2coef' t1 v)))); neuper@37950: t2pp:=(#2(hd(the(term2coef' t2 v)))); neuper@37950: t1c:=(#1(hd(the(term2coef' t1 v)))); neuper@37950: t2c:=(#1(hd(the(term2coef' t2 v)))); neuper@37950: neuper@37950: SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] handle _ => NONE neuper@37950: neuper@37950: ) neuper@37950: end neuper@37950: | term2coef' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ (t2 as Free (str2,_))) v :mv_poly option= neuper@37950: let neuper@38006: val x= Unsynchronized.ref NONE; neuper@38006: val len= Unsynchronized.ref 0; neuper@38006: val vl= Unsynchronized.ref []; neuper@38006: val vh= Unsynchronized.ref []; neuper@38006: val vtemp= Unsynchronized.ref []; neuper@38006: val i= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: if (not o is_numeral) str1 andalso is_numeral str2 then neuper@37950: ( neuper@37950: len:=length(v); neuper@37950: vh:=v; neuper@37950: neuper@37950: while ((!len)>(!i)) do neuper@37950: ( neuper@37950: if str1=hd((!vh)) then neuper@37950: ( neuper@37950: vl:=((the o int_of_str) str2)::(!vl) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: vl:=0::(!vl) neuper@37950: ); neuper@37950: vh:=tl(!vh); neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: SOME [(1,rev(!vl))] handle _ => NONE neuper@37950: ) neuper@37950: else raise error ("RATIONALS_TERM2COEF_EXCEPTION 1: Invalid term") neuper@37950: ) neuper@37950: end neuper@38014: | term2coef' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option= neuper@37950: ( neuper@37950: SOME ((the(term2coef' t1 v)) @ (the(term2coef' t2 v))) handle _ => NONE neuper@37950: ) neuper@38014: | term2coef' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option= neuper@37950: ( neuper@37950: SOME ((the(term2coef' t1 v)) @ mv_skalar_mul((the(term2coef' t2 v)),1)) handle _ => NONE neuper@37950: ) neuper@37950: | term2coef' (term) v = raise error ("RATIONALS_TERM2COEF_EXCEPTION 2: Invalid term"); neuper@37950: neuper@37950: (*. checks if all coefficients of a polynomial are positiv (except the first) .*) neuper@37950: fun check_coeff t = (* erste Koeffizient kann <0 sein !!! *) neuper@37950: if count_neg(tl(the(term2coef' t (get_vars(t)))))=0 then true neuper@37950: else false; neuper@37950: neuper@37950: (*. checks for expanded term [3] .*) neuper@37950: fun is_expanded t = test_exp t " " andalso check_coeff(t); neuper@37950: neuper@37950: (*WN.7.3.03 Hilfsfunktion f"ur term2poly'*) neuper@37950: fun mk_monom v' p vs = neuper@37950: let fun conv p (v: string) = if v'= v then p else 0 neuper@37950: in map (conv p) vs end; neuper@37950: (* mk_monom "y" 5 ["a","b","x","y","z"]; neuper@37950: val it = [0,0,0,5,0] : int list*) neuper@37950: neuper@37950: (*. this function converts the term representation into the internal representation mv_poly .*) neuper@37950: fun term2poly' (Const ("uminus",_) $ Free (str,_)) v = (*WN.7.3.03*) neuper@37950: if is_numeral str neuper@37950: then SOME [((the o int_of_str) ("-"^str), mk_monom "#" 0 v)] neuper@37950: else SOME [(~1, mk_monom str 1 v)] neuper@37950: neuper@37950: | term2poly' (Free(str,_)) v :mv_poly option = neuper@37950: let neuper@38006: val x= Unsynchronized.ref NONE; neuper@38006: val len= Unsynchronized.ref 0; neuper@38006: val vl= Unsynchronized.ref []; neuper@38006: val vh= Unsynchronized.ref []; neuper@38006: val i= Unsynchronized.ref 0; neuper@37950: in neuper@37950: if is_numeral str then neuper@37950: ( neuper@37950: SOME [(((the o int_of_str) str),mv_null2 v)] handle _ => NONE neuper@37950: ) neuper@37950: else (* variable *) neuper@37950: ( neuper@37950: len:=length v; neuper@37950: vh:= v; neuper@37950: while ((!len)>(!i)) do neuper@37950: ( neuper@37950: if str=hd((!vh)) then neuper@37950: ( neuper@37950: vl:=1::(!vl) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: vl:=0::(!vl) neuper@37950: ); neuper@37950: vh:=tl(!vh); neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: SOME [(1,rev(!vl))] handle _ => NONE neuper@37950: ) neuper@37950: end neuper@37950: | term2poly' (Const ("op *",_) $ t1 $ t2) v :mv_poly option= neuper@37950: let neuper@38006: val t1pp= Unsynchronized.ref []; neuper@38006: val t2pp= Unsynchronized.ref []; neuper@38006: val t1c= Unsynchronized.ref 0; neuper@38006: val t2c= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: t1pp:=(#2(hd(the(term2poly' t1 v)))); neuper@37950: t2pp:=(#2(hd(the(term2poly' t2 v)))); neuper@37950: t1c:=(#1(hd(the(term2poly' t1 v)))); neuper@37950: t2c:=(#1(hd(the(term2poly' t2 v)))); neuper@37950: neuper@37950: SOME [( (!t1c)*(!t2c) ,( (map op+) ((!t1pp)~~(!t2pp)) ) )] neuper@37950: handle _ => NONE neuper@37950: neuper@37950: ) neuper@37950: end neuper@37950: | term2poly' (Const ("Atools.pow",_) $ (t1 as Free (str1,_)) $ neuper@37950: (t2 as Free (str2,_))) v :mv_poly option= neuper@37950: let neuper@38006: val x= Unsynchronized.ref NONE; neuper@38006: val len= Unsynchronized.ref 0; neuper@38006: val vl= Unsynchronized.ref []; neuper@38006: val vh= Unsynchronized.ref []; neuper@38006: val vtemp= Unsynchronized.ref []; neuper@38006: val i= Unsynchronized.ref 0; neuper@37950: in neuper@37950: ( neuper@37950: if (not o is_numeral) str1 andalso is_numeral str2 then neuper@37950: ( neuper@37950: len:=length(v); neuper@37950: vh:=v; neuper@37950: neuper@37950: while ((!len)>(!i)) do neuper@37950: ( neuper@37950: if str1=hd((!vh)) then neuper@37950: ( neuper@37950: vl:=((the o int_of_str) str2)::(!vl) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: vl:=0::(!vl) neuper@37950: ); neuper@37950: vh:=tl(!vh); neuper@37950: i:=(!i)+1 neuper@37950: ); neuper@37950: SOME [(1,rev(!vl))] handle _ => NONE neuper@37950: ) neuper@37950: else raise error ("RATIONALS_TERM2POLY_EXCEPTION 1: Invalid term") neuper@37950: ) neuper@37950: end neuper@38014: | term2poly' (Const ("Groups.plus_class.plus",_) $ t1 $ t2) v :mv_poly option = neuper@37950: ( neuper@37950: SOME ((the(term2poly' t1 v)) @ (the(term2poly' t2 v))) handle _ => NONE neuper@37950: ) neuper@38014: | term2poly' (Const ("Groups.minus_class.minus",_) $ t1 $ t2) v :mv_poly option = neuper@37950: ( neuper@37950: SOME ((the(term2poly' t1 v)) @ mv_skalar_mul((the(term2poly' t2 v)),~1)) handle _ => NONE neuper@37950: ) neuper@37950: | term2poly' (term) v = raise error ("RATIONALS_TERM2POLY_EXCEPTION 2: Invalid term"); neuper@37950: neuper@37950: (*. translates an Isabelle term into internal representation. neuper@37950: term2poly neuper@37950: fn : term -> (*normalform [2] *) neuper@37950: string list -> (*for ...!!! BITTE DIE ERKLÄRUNG, neuper@37950: DIE DU MIR LETZTES MAL GEGEBEN HAST*) neuper@37950: mv_monom list (*internal representation *) neuper@37950: option (*the translation may fail with NONE*) neuper@37950: .*) neuper@37950: fun term2poly (t:term) v = neuper@37950: if is_polynomial t then term2poly' t v neuper@37950: else raise error ("term2poly: invalid = "^(term2str t)); neuper@37950: neuper@37950: (*. same as term2poly with automatic detection of the variables .*) neuper@37950: fun term2polyx t = term2poly t (((map free2str) o vars) t); neuper@37950: neuper@37950: (*. checks if the term is in expanded polynomial form and converts it into the internal representation .*) neuper@37950: fun expanded2poly (t:term) v = neuper@37950: (*if is_expanded t then*) term2poly' t v neuper@37950: (*else raise error ("RATIONALS_EXPANDED2POLY_EXCEPTION: Invalid Polynomial")*); neuper@37950: neuper@37950: (*. same as expanded2poly with automatic detection of the variables .*) neuper@37950: fun expanded2polyx t = expanded2poly t (((map free2str) o vars) t); neuper@37950: neuper@37950: (*. converts a powerproduct into term representation .*) neuper@37950: fun powerproduct2term(xs,v) = neuper@37950: let neuper@38006: val xss= Unsynchronized.ref xs; neuper@38006: val vv= Unsynchronized.ref v; neuper@37950: in neuper@37950: ( neuper@37950: while hd(!xss)=0 do neuper@37950: ( neuper@37950: xss:=tl(!xss); neuper@37950: vv:=tl(!vv) neuper@37950: ); neuper@37950: neuper@37950: if list_is_null(tl(!xss)) then neuper@37950: ( neuper@37950: if hd(!xss)=1 then Free(hd(!vv), HOLogic.realT) neuper@37950: else neuper@37950: ( neuper@37950: Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT) neuper@37950: ) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if hd(!xss)=1 then neuper@37950: ( neuper@37950: Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: Free(hd(!vv), HOLogic.realT) $ neuper@37950: powerproduct2term(tl(!xss),tl(!vv)) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: ( neuper@37950: Const("Atools.pow",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: Free(hd(!vv), HOLogic.realT) $ Free(str_of_int (hd(!xss)),HOLogic.realT) neuper@37950: ) $ neuper@37950: powerproduct2term(tl(!xss),tl(!vv)) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. converts a monom into term representation .*) neuper@37950: (*fun monom2term ((c,e):mv_monom, v:string list) = neuper@37950: if c=0 then Free(str_of_int 0,HOLogic.realT) neuper@37950: else neuper@37950: ( neuper@37950: if list_is_null(e) then neuper@37950: ( neuper@37950: Free(str_of_int c,HOLogic.realT) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if c=1 then neuper@37950: ( neuper@37950: powerproduct2term(e,v) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: Free(str_of_int c,HOLogic.realT) $ neuper@37950: powerproduct2term(e,v) neuper@37950: ) neuper@37950: ) neuper@37950: );*) neuper@37950: neuper@37950: neuper@37950: (*fun monom2term ((i, is):mv_monom, v) = neuper@37950: if list_is_null is neuper@37950: then neuper@37950: if i >= 0 neuper@37950: then Free (str_of_int i, HOLogic.realT) neuper@37950: else Const ("uminus", HOLogic.realT --> HOLogic.realT) $ neuper@37950: Free ((str_of_int o abs) i, HOLogic.realT) neuper@37950: else neuper@37950: if i > 0 neuper@37950: then Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $ neuper@37950: (Free (str_of_int i, HOLogic.realT)) $ neuper@37950: powerproduct2term(is, v) neuper@37950: else Const ("op *", [HOLogic.realT,HOLogic.realT]---> HOLogic.realT) $ neuper@37950: (Const ("uminus", HOLogic.realT --> HOLogic.realT) $ neuper@37950: Free ((str_of_int o abs) i, HOLogic.realT)) $ neuper@37950: powerproduct2term(is, vs);---------------------------*) neuper@37950: fun monom2term ((i, is) : mv_monom, vs) = neuper@37950: if list_is_null is neuper@37950: then Free (str_of_int i, HOLogic.realT) neuper@37950: else if i = 1 neuper@37950: then powerproduct2term (is, vs) neuper@37950: else Const ("op *", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $ neuper@37950: (Free (str_of_int i, HOLogic.realT)) $ neuper@37950: powerproduct2term (is, vs); neuper@37950: neuper@37950: (*. converts the internal polynomial representation into an Isabelle term.*) neuper@37950: fun poly2term' ([] : mv_poly, vs) = Free(str_of_int 0, HOLogic.realT) neuper@37950: | poly2term' ([(c, e) : mv_monom], vs) = monom2term ((c, e), vs) neuper@37950: | poly2term' ((c, e) :: ces, vs) = neuper@38014: Const("Groups.plus_class.plus", [HOLogic.realT, HOLogic.realT] ---> HOLogic.realT) $ neuper@37950: poly2term (ces, vs) $ monom2term ((c, e), vs) neuper@37950: and poly2term (xs, vs) = poly2term' (rev (sort (mv_geq LEX_) (xs)), vs); neuper@37950: neuper@37950: neuper@37950: (*. converts a monom into term representation .*) neuper@37950: (*. ignores the sign of the coefficients => use only for exp-poly functions .*) neuper@37950: fun monom2term2((c,e):mv_monom, v:string list) = neuper@37950: if c=0 then Free(str_of_int 0,HOLogic.realT) neuper@37950: else neuper@37950: ( neuper@37950: if list_is_null(e) then neuper@37950: ( neuper@37950: Free(str_of_int (abs(c)),HOLogic.realT) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: if abs(c)=1 then neuper@37950: ( neuper@37950: powerproduct2term(e,v) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: Const("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: Free(str_of_int (abs(c)),HOLogic.realT) $ neuper@37950: powerproduct2term(e,v) neuper@37950: ) neuper@37950: ) neuper@37950: ); neuper@37950: neuper@37950: (*. converts the expanded polynomial representation into the term representation .*) neuper@37950: fun exp2term' ([]:mv_poly,vars) = Free(str_of_int 0,HOLogic.realT) neuper@37950: | exp2term' ([(c,e)],vars) = monom2term((c,e),vars) neuper@37950: | exp2term' ((c1,e1)::others,vars) = neuper@37950: if c1<0 then neuper@38014: Const("Groups.minus_class.minus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: exp2term'(others,vars) $ neuper@37950: ( neuper@37950: monom2term2((c1,e1),vars) neuper@37950: ) neuper@37950: else neuper@38014: Const("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: exp2term'(others,vars) $ neuper@37950: ( neuper@37950: monom2term2((c1,e1),vars) neuper@37950: ); neuper@37950: neuper@37950: (*. sorts the powerproduct by lexicographic termorder and converts them into neuper@37950: a term in polynomial representation .*) neuper@37950: fun poly2expanded (xs,vars) = exp2term'(rev(sort (mv_geq LEX_) (xs)),vars); neuper@37950: neuper@37950: (*. converts a polynomial into expanded form .*) neuper@37950: fun polynomial2expanded t = neuper@37950: (let neuper@37950: val vars=(((map free2str) o vars) t); neuper@37950: in neuper@37950: SOME (poly2expanded (the (term2poly t vars), vars)) neuper@37950: end) handle _ => NONE; neuper@37950: neuper@37950: (*. converts a polynomial into polynomial form .*) neuper@37950: fun expanded2polynomial t = neuper@37950: (let neuper@37950: val vars=(((map free2str) o vars) t); neuper@37950: in neuper@37950: SOME (poly2term (the (expanded2poly t vars), vars)) neuper@37950: end) handle _ => NONE; neuper@37950: neuper@37950: neuper@37950: (*. calculates the greatest common divisor of numerator and denominator and seperates it from each .*) neuper@38014: fun step_cancel (t as Const ("Rings.inverse_class.divide",_) $ p1 $ p2) = neuper@37950: let neuper@38006: val p1' = Unsynchronized.ref []; neuper@38006: val p2' = Unsynchronized.ref []; neuper@38006: val p3 = Unsynchronized.ref [] neuper@37950: val vars = rev(get_vars(p1) union get_vars(p2)); neuper@37950: in neuper@37950: ( neuper@37950: p1':= sort (mv_geq LEX_) (the (term2poly p1 vars )); neuper@37950: p2':= sort (mv_geq LEX_) (the (term2poly p2 vars )); neuper@37950: p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2')); neuper@37950: if (!p3)=[(1,mv_null2(vars))] then neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2 neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: neuper@37950: p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_))); neuper@37950: p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_))); neuper@37950: neuper@37950: if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(!p1',vars) $ neuper@37950: poly2term(!p3,vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(!p2',vars) $ neuper@37950: poly2term(!p3,vars) neuper@37950: ) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=mv_skalar_mul(!p1',~1); neuper@37950: p2':=mv_skalar_mul(!p2',~1); neuper@37950: p3:=mv_skalar_mul(!p3,~1); neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(!p1',vars) $ neuper@37950: poly2term(!p3,vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(!p2',vars) $ neuper@37950: poly2term(!p3,vars) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | step_cancel _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction"); neuper@37950: neuper@37950: neuper@37950: (*. same as step_cancel, this time for expanded forms (input+output) .*) neuper@38014: fun step_cancel_expanded (t as Const ("Rings.inverse_class.divide",_) $ p1 $ p2) = neuper@37950: let neuper@38006: val p1' = Unsynchronized.ref []; neuper@38006: val p2' = Unsynchronized.ref []; neuper@38006: val p3 = Unsynchronized.ref [] neuper@37950: val vars = rev(get_vars(p1) union get_vars(p2)); neuper@37950: in neuper@37950: ( neuper@37950: p1':= sort (mv_geq LEX_) (the (expanded2poly p1 vars )); neuper@37950: p2':= sort (mv_geq LEX_) (the (expanded2poly p2 vars )); neuper@37950: p3:= sort (mv_geq LEX_) (mv_gcd (!p1') (!p2')); neuper@37950: if (!p3)=[(1,mv_null2(vars))] then neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2 neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: neuper@37950: p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_))); neuper@37950: p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_))); neuper@37950: neuper@37950: if #1(hd(sort (mv_geq LEX_) (!p2')))(* mv_lc2(!p2',LEX_)*)>0 then neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(!p1',vars) $ neuper@37950: poly2expanded(!p3,vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(!p2',vars) $ neuper@37950: poly2expanded(!p3,vars) neuper@37950: ) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=mv_skalar_mul(!p1',~1); neuper@37950: p2':=mv_skalar_mul(!p2',~1); neuper@37950: p3:=mv_skalar_mul(!p3,~1); neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(!p1',vars) $ neuper@37950: poly2expanded(!p3,vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(!p2',vars) $ neuper@37950: poly2expanded(!p3,vars) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | step_cancel_expanded _ = raise error ("RATIONALS_STEP_CANCEL_EXCEPTION: Invalid fraction"); neuper@37950: neuper@37950: (*. calculates the greatest common divisor of numerator and denominator and divides each through it .*) neuper@38014: fun direct_cancel (t as Const ("Rings.inverse_class.divide",_) $ p1 $ p2) = neuper@37950: let neuper@38006: val p1' = Unsynchronized.ref []; neuper@38006: val p2' = Unsynchronized.ref []; neuper@38006: val p3 = Unsynchronized.ref [] neuper@37950: val vars = rev(get_vars(p1) union get_vars(p2)); neuper@37950: in neuper@37950: ( neuper@37950: p1':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p1 vars )),LEX_)); neuper@37950: p2':=sort (mv_geq LEX_) (mv_shorten((the (term2poly p2 vars )),LEX_)); neuper@37950: p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2')); neuper@37950: neuper@37950: if (!p3)=[(1,mv_null2(vars))] then neuper@37950: ( neuper@38014: (Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[]) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_))); neuper@37950: p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_))); neuper@37950: if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then neuper@37950: ( neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!p1'),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!p2'),vars) neuper@37950: ) neuper@37950: ) neuper@37950: , neuper@37950: if mv_grad(!p3)>0 then neuper@37950: [ neuper@37950: ( neuper@37950: Const ("Not",[bool]--->bool) $ neuper@37950: ( neuper@37950: Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $ neuper@37950: poly2term((!p3),vars) $ neuper@37950: Free("0",HOLogic.realT) neuper@37950: ) neuper@37950: ) neuper@37950: ] neuper@37950: else neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=mv_skalar_mul(!p1',~1); neuper@37950: p2':=mv_skalar_mul(!p2',~1); neuper@37950: if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1); neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!p1'),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!p2'),vars) neuper@37950: ) neuper@37950: , neuper@37950: if mv_grad(!p3)>0 then neuper@37950: [ neuper@37950: ( neuper@37950: Const ("Not",[bool]--->bool) $ neuper@37950: ( neuper@37950: Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $ neuper@37950: poly2term((!p3),vars) $ neuper@37950: Free("0",HOLogic.realT) neuper@37950: ) neuper@37950: ) neuper@37950: ] neuper@37950: else neuper@37950: [] neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | direct_cancel _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction"); neuper@37950: neuper@37950: (*. same es direct_cancel, this time for expanded forms (input+output).*) neuper@38014: fun direct_cancel_expanded (t as Const ("Rings.inverse_class.divide",_) $ p1 $ p2) = neuper@37950: let neuper@38006: val p1' = Unsynchronized.ref []; neuper@38006: val p2' = Unsynchronized.ref []; neuper@38006: val p3 = Unsynchronized.ref [] neuper@37950: val vars = rev(get_vars(p1) union get_vars(p2)); neuper@37950: in neuper@37950: ( neuper@37950: p1':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p1 vars )),LEX_)); neuper@37950: p2':=sort (mv_geq LEX_) (mv_shorten((the (expanded2poly p2 vars )),LEX_)); neuper@37950: p3 :=sort (mv_geq LEX_) (mv_gcd (!p1') (!p2')); neuper@37950: neuper@37950: if (!p3)=[(1,mv_null2(vars))] then neuper@37950: ( neuper@38014: (Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ p1 $ p2,[]) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=sort (mv_geq LEX_) (#1(mv_division((!p1'),(!p3),LEX_))); neuper@37950: p2':=sort (mv_geq LEX_) (#1(mv_division((!p2'),(!p3),LEX_))); neuper@37950: if #1(hd(sort (mv_geq LEX_) (!p2'))) (*mv_lc2(!p2',LEX_)*)>0 then neuper@37950: ( neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!p1'),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!p2'),vars) neuper@37950: ) neuper@37950: ) neuper@37950: , neuper@37950: if mv_grad(!p3)>0 then neuper@37950: [ neuper@37950: ( neuper@37950: Const ("Not",[bool]--->bool) $ neuper@37950: ( neuper@37950: Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $ neuper@37950: poly2expanded((!p3),vars) $ neuper@37950: Free("0",HOLogic.realT) neuper@37950: ) neuper@37950: ) neuper@37950: ] neuper@37950: else neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: p1':=mv_skalar_mul(!p1',~1); neuper@37950: p2':=mv_skalar_mul(!p2',~1); neuper@37950: if length(!p3)> 2*(count_neg(!p3)) then () else p3 :=mv_skalar_mul(!p3,~1); neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!p1'),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!p2'),vars) neuper@37950: ) neuper@37950: , neuper@37950: if mv_grad(!p3)>0 then neuper@37950: [ neuper@37950: ( neuper@37950: Const ("Not",[bool]--->bool) $ neuper@37950: ( neuper@37950: Const("op =",[HOLogic.realT,HOLogic.realT]--->bool) $ neuper@37950: poly2expanded((!p3),vars) $ neuper@37950: Free("0",HOLogic.realT) neuper@37950: ) neuper@37950: ) neuper@37950: ] neuper@37950: else neuper@37950: [] neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | direct_cancel_expanded _ = raise error ("RATIONALS_DIRECT_CANCEL_EXCEPTION: Invalid fraction"); neuper@37950: neuper@37950: neuper@37950: (*. adds two fractions .*) neuper@38014: fun add_fract ((Const("Rings.inverse_class.divide",_) $ t11 $ t12),(Const("Rings.inverse_class.divide",_) $ t21 $ t22)) = neuper@37950: let neuper@37950: val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22); neuper@38006: val t11'= Unsynchronized.ref (the(term2poly t11 vars)); neuper@37950: val _= writeln"### add_fract: done t11" neuper@38006: val t12'= Unsynchronized.ref (the(term2poly t12 vars)); neuper@37950: val _= writeln"### add_fract: done t12" neuper@38006: val t21'= Unsynchronized.ref (the(term2poly t21 vars)); neuper@37950: val _= writeln"### add_fract: done t21" neuper@38006: val t22'= Unsynchronized.ref (the(term2poly t22 vars)); neuper@37950: val _= writeln"### add_fract: done t22" neuper@38006: val den= Unsynchronized.ref []; neuper@38006: val nom= Unsynchronized.ref []; neuper@38006: val m1= Unsynchronized.ref []; neuper@38006: val m2= Unsynchronized.ref []; neuper@37950: in neuper@37950: neuper@37950: ( neuper@37950: den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22')); neuper@37950: writeln"### add_fract: done sort mv_lcm"; neuper@37950: m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_))); neuper@37950: writeln"### add_fract: done sort mv_division t12"; neuper@37950: m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_))); neuper@37950: writeln"### add_fract: done sort mv_division t22"; neuper@37950: nom :=sort (mv_geq LEX_) neuper@37950: (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_), neuper@37950: mv_mul(!t21',!m2,LEX_), neuper@37950: LEX_), neuper@37950: LEX_)); neuper@37950: writeln"### add_fract: done sort mv_add"; neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!nom),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2term((!den),vars) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | add_fract (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: Invalid add_fraction call"); neuper@37950: neuper@37950: (*. adds two expanded fractions .*) neuper@38014: fun add_fract_exp ((Const("Rings.inverse_class.divide",_) $ t11 $ t12),(Const("Rings.inverse_class.divide",_) $ t21 $ t22)) = neuper@37950: let neuper@37950: val vars=get_vars(t11) union get_vars(t12) union get_vars(t21) union get_vars(t22); neuper@38006: val t11'= Unsynchronized.ref (the(expanded2poly t11 vars)); neuper@38006: val t12'= Unsynchronized.ref (the(expanded2poly t12 vars)); neuper@38006: val t21'= Unsynchronized.ref (the(expanded2poly t21 vars)); neuper@38006: val t22'= Unsynchronized.ref (the(expanded2poly t22 vars)); neuper@38006: val den= Unsynchronized.ref []; neuper@38006: val nom= Unsynchronized.ref []; neuper@38006: val m1= Unsynchronized.ref []; neuper@38006: val m2= Unsynchronized.ref []; neuper@37950: in neuper@37950: neuper@37950: ( neuper@37950: den :=sort (mv_geq LEX_) (mv_lcm (!t12') (!t22')); neuper@37950: m1 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t12',LEX_))); neuper@37950: m2 :=sort (mv_geq LEX_) (#1(mv_division(!den,!t22',LEX_))); neuper@37950: nom :=sort (mv_geq LEX_) (mv_shorten(mv_add(mv_mul(!t11',!m1,LEX_),mv_mul(!t21',!m2,LEX_),LEX_),LEX_)); neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!nom),vars) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: poly2expanded((!den),vars) neuper@37950: ) neuper@37950: ) neuper@37950: ) neuper@37950: end neuper@37950: | add_fract_exp (_,_) = raise error ("RATIONALS_ADD_FRACTION_EXP_EXCEPTION: Invalid add_fraction call"); neuper@37950: neuper@37950: (*. adds a list of terms .*) neuper@37950: fun add_list_of_fractions []= (Free("0",HOLogic.realT),[]) neuper@37950: | add_list_of_fractions [x]= direct_cancel x neuper@37950: | add_list_of_fractions (x::y::xs) = neuper@37950: let neuper@37950: val (t1a,rest1)=direct_cancel(x); neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(x)"; neuper@37950: val (t2a,rest2)=direct_cancel(y); neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(y)"; neuper@37950: val (t3a,rest3)=(add_list_of_fractions (add_fract(t1a,t2a)::xs)); neuper@37950: val _= writeln"### add_list_of_fractions xs: has done add_list_of_fraction xs"; neuper@37950: val (t4a,rest4)=direct_cancel(t3a); neuper@37950: val _= writeln"### add_list_of_fractions xs: has done direct_cancel(t3a)"; neuper@37950: val rest=rest1 union rest2 union rest3 union rest4; neuper@37950: in neuper@37950: (writeln"### add_list_of_fractions in"; neuper@37950: ( neuper@37950: (t4a,rest) neuper@37950: ) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. adds a list of expanded terms .*) neuper@37950: fun add_list_of_fractions_exp []= (Free("0",HOLogic.realT),[]) neuper@37950: | add_list_of_fractions_exp [x]= direct_cancel_expanded x neuper@37950: | add_list_of_fractions_exp (x::y::xs) = neuper@37950: let neuper@37950: val (t1a,rest1)=direct_cancel_expanded(x); neuper@37950: val (t2a,rest2)=direct_cancel_expanded(y); neuper@37950: val (t3a,rest3)=(add_list_of_fractions_exp (add_fract_exp(t1a,t2a)::xs)); neuper@37950: val (t4a,rest4)=direct_cancel_expanded(t3a); neuper@37950: val rest=rest1 union rest2 union rest3 union rest4; neuper@37950: in neuper@37950: ( neuper@37950: (t4a,rest) neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the lcm of a list of mv_poly .*) neuper@37950: fun calc_lcm ([x],var)= (x,var) neuper@37950: | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var); neuper@37950: neuper@37950: (*. converts a list of terms to a list of mv_poly .*) neuper@37950: fun t2d([],_)=[] neuper@38014: | t2d((t as (Const("Rings.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars); neuper@37950: neuper@37950: (*. same as t2d, this time for expanded forms .*) neuper@37950: fun t2d_exp([],_)=[] neuper@38014: | t2d_exp((t as (Const("Rings.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars); neuper@37950: neuper@37950: (*. converts a list of fract terms to a list of their denominators .*) neuper@37950: fun termlist2denominators [] = ([],[]) neuper@37950: | termlist2denominators xs = neuper@37950: let neuper@38006: val xxs= Unsynchronized.ref xs; neuper@38006: val var= Unsynchronized.ref []; neuper@37950: in neuper@37950: var:=[]; neuper@37950: while length(!xxs)>0 do neuper@37950: ( neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs); neuper@37950: in neuper@37950: ( neuper@37950: xxs:=tl(!xxs); neuper@37950: var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var)) neuper@37950: ) neuper@37950: end neuper@37950: ); neuper@37950: (t2d(xs,!var),!var) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the lcm of a list of mv_poly .*) neuper@37950: fun calc_lcm ([x],var)= (x,var) neuper@37950: | calc_lcm ((x::xs),var) = (mv_lcm x (#1(calc_lcm (xs,var))),var); neuper@37950: neuper@37950: (*. converts a list of terms to a list of mv_poly .*) neuper@37950: fun t2d([],_)=[] neuper@38014: | t2d((t as (Const("Rings.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(term2poly p2 vars)) :: t2d(xs,vars); neuper@37950: neuper@37950: (*. same as t2d, this time for expanded forms .*) neuper@37950: fun t2d_exp([],_)=[] neuper@38014: | t2d_exp((t as (Const("Rings.inverse_class.divide",_) $ p1 $ p2))::xs,vars)= (the(expanded2poly p2 vars)) :: t2d_exp(xs,vars); neuper@37950: neuper@37950: (*. converts a list of fract terms to a list of their denominators .*) neuper@37950: fun termlist2denominators [] = ([],[]) neuper@37950: | termlist2denominators xs = neuper@37950: let neuper@38006: val xxs= Unsynchronized.ref xs; neuper@38006: val var= Unsynchronized.ref []; neuper@37950: in neuper@37950: var:=[]; neuper@37950: while length(!xxs)>0 do neuper@37950: ( neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs); neuper@37950: in neuper@37950: ( neuper@37950: xxs:=tl(!xxs); neuper@37950: var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var)) neuper@37950: ) neuper@37950: end neuper@37950: ); neuper@37950: (t2d(xs,!var),!var) neuper@37950: end; neuper@37950: neuper@37950: (*. same as termlist2denminators, this time for expanded forms .*) neuper@37950: fun termlist2denominators_exp [] = ([],[]) neuper@37950: | termlist2denominators_exp xs = neuper@37950: let neuper@38006: val xxs= Unsynchronized.ref xs; neuper@38006: val var= Unsynchronized.ref []; neuper@37950: in neuper@37950: var:=[]; neuper@37950: while length(!xxs)>0 do neuper@37950: ( neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1x $ p2x)=hd(!xxs); neuper@37950: in neuper@37950: ( neuper@37950: xxs:=tl(!xxs); neuper@37950: var:=((get_vars(p2x)) union (get_vars(p1x)) union (!var)) neuper@37950: ) neuper@37950: end neuper@37950: ); neuper@37950: (t2d_exp(xs,!var),!var) neuper@37950: end; neuper@37950: neuper@37950: (*. reduces all fractions to the least common denominator .*) neuper@37950: fun com_den(x::xs,denom,den,var)= neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1' $ p2')=x; neuper@37950: val p2= sort (mv_geq LEX_) (the(term2poly p2' var)); neuper@37950: val p3= #1(mv_division(denom,p2,LEX_)); neuper@37950: val p1var=get_vars(p1'); neuper@37950: in neuper@37950: if length(xs)>0 then neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: poly2term(the (term2poly p1' p1var),p1var) neuper@37950: $ neuper@37950: den neuper@37950: ) neuper@37950: $ neuper@37950: #1(com_den(xs,denom,den,var)) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(the (term2poly p1' p1var),p1var) $ neuper@37950: poly2term(p3,var) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: den neuper@37950: ) neuper@37950: ) neuper@37950: $ neuper@37950: #1(com_den(xs,denom,den,var)) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: poly2term(the (term2poly p1' p1var),p1var) neuper@37950: $ neuper@37950: den neuper@37950: ) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(the (term2poly p1' p1var),p1var) $ neuper@37950: poly2term(p3,var) neuper@37950: ) neuper@37950: $ neuper@37950: den neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (*. same as com_den, this time for expanded forms .*) neuper@37950: fun com_den_exp(x::xs,denom,den,var)= neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1' $ p2')=x; neuper@37950: val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var)); neuper@37950: val p3= #1(mv_division(denom,p2,LEX_)); neuper@37950: val p1var=get_vars(p1'); neuper@37950: in neuper@37950: if length(xs)>0 then neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: poly2expanded(the(expanded2poly p1' p1var),p1var) neuper@37950: $ neuper@37950: den neuper@37950: ) neuper@37950: $ neuper@37950: #1(com_den_exp(xs,denom,den,var)) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(the(expanded2poly p1' p1var),p1var) $ neuper@37950: poly2expanded(p3,var) neuper@37950: ) neuper@37950: $ neuper@37950: ( neuper@37950: den neuper@37950: ) neuper@37950: ) neuper@37950: $ neuper@37950: #1(com_den_exp(xs,denom,den,var)) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: poly2expanded(the(expanded2poly p1' p1var),p1var) neuper@37950: $ neuper@37950: den neuper@37950: ) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) neuper@37950: $ neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(the(expanded2poly p1' p1var),p1var) $ neuper@37950: poly2expanded(p3,var) neuper@37950: ) neuper@37950: $ neuper@37950: den neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (* wird aktuell nicht mehr gebraucht, bei rückänderung schon neuper@37950: ------------------------------------------------------------- neuper@37950: (* WN0210???SK brauch ma des überhaupt *) neuper@37950: fun com_den2(x::xs,denom,den,var)= neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1' $ p2')=x; neuper@37950: val p2= sort (mv_geq LEX_) (the(term2poly p2' var)); neuper@37950: val p3= #1(mv_division(denom,p2,LEX_)); neuper@37950: val p1var=get_vars(p1'); neuper@37950: in neuper@37950: if length(xs)>0 then neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(the(term2poly p1' p1var),p1var) $ neuper@37950: com_den2(xs,denom,den,var) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: ( neuper@37950: let neuper@37950: val p3'=poly2term(p3,var); neuper@37950: val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3'); neuper@37950: in neuper@37950: poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars) neuper@37950: end neuper@37950: ) $ neuper@37950: com_den2(xs,denom,den,var) neuper@37950: ) neuper@37950: else neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@37950: poly2term(the(term2poly p1' p1var),p1var) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: let neuper@37950: val p3'=poly2term(p3,var); neuper@37950: val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3'); neuper@37950: in neuper@37950: poly2term(sort (mv_geq LEX_) (mv_mul(the(term2poly p1' vars) ,the(term2poly p3' vars),LEX_)),vars) neuper@37950: end neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (* WN0210???SK brauch ma des überhaupt *) neuper@37950: fun com_den_exp2(x::xs,denom,den,var)= neuper@37950: let neuper@38014: val (t as Const ("Rings.inverse_class.divide",_) $ p1' $ p2')=x; neuper@37950: val p2= sort (mv_geq LEX_) (the(expanded2poly p2' var)); neuper@37950: val p3= #1(mv_division(denom,p2,LEX_)); neuper@37950: val p1var=get_vars p1'; neuper@37950: in neuper@37950: if length(xs)>0 then neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(the (expanded2poly p1' p1var),p1var) $ neuper@37950: com_den_exp2(xs,denom,den,var) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@38014: Const ("Groups.plus_class.plus",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: ( neuper@37950: let neuper@37950: val p3'=poly2expanded(p3,var); neuper@37950: val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3'); neuper@37950: in neuper@37950: poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars) neuper@37950: end neuper@37950: ) $ neuper@37950: com_den_exp2(xs,denom,den,var) neuper@37950: ) neuper@37950: else neuper@37950: if p3=[(1,mv_null2(var))] then neuper@37950: ( neuper@37950: poly2expanded(the (expanded2poly p1' p1var),p1var) neuper@37950: ) neuper@37950: else neuper@37950: ( neuper@37950: let neuper@37950: val p3'=poly2expanded(p3,var); neuper@37950: val vars= (((map free2str) o vars) p1') union (((map free2str) o vars) p3'); neuper@37950: in neuper@37950: poly2expanded(sort (mv_geq LEX_) (mv_mul(the(expanded2poly p1' vars) ,the(expanded2poly p3' vars),LEX_)),vars) neuper@37950: end neuper@37950: ) neuper@37950: end; neuper@37950: ---------------------------------------------------------*) neuper@37950: neuper@37950: neuper@37950: (*. searches for an element y of a list ys, which has an gcd not 1 with x .*) neuper@37950: fun exists_gcd (x,[]) = false neuper@37950: | exists_gcd (x,y::ys) = if mv_gcd x y = [(1,mv_null2(#2(hd(x))))] then exists_gcd (x,ys) neuper@37950: else true; neuper@37950: neuper@37950: (*. divides each element of the list xs with y .*) neuper@37950: fun list_div ([],y) = [] neuper@37950: | list_div (x::xs,y) = neuper@37950: let neuper@37950: val (d,r)=mv_division(x,y,LEX_); neuper@37950: in neuper@37950: if r=[] then neuper@37950: d::list_div(xs,y) neuper@37950: else x::list_div(xs,y) neuper@37950: end; neuper@37950: neuper@37950: (*. checks if x is in the list ys .*) neuper@37950: fun in_list (x,[]) = false neuper@37950: | in_list (x,y::ys) = if x=y then true neuper@37950: else in_list(x,ys); neuper@37950: neuper@37950: (*. deletes all equal elements of the list xs .*) neuper@37950: fun kill_equal [] = [] neuper@37950: | kill_equal (x::xs) = if in_list(x,xs) orelse x=[(1,mv_null2(#2(hd(x))))] then kill_equal(xs) neuper@37950: else x::kill_equal(xs); neuper@37950: neuper@37950: (*. searches for new factors .*) neuper@37950: fun new_factors [] = [] neuper@37950: | new_factors (list:mv_poly list):mv_poly list = neuper@37950: let neuper@37950: val l = kill_equal list; neuper@37950: val len = length(l); neuper@37950: in neuper@37950: if len>=2 then neuper@37950: ( neuper@37950: let neuper@37950: val x::y::xs=l; neuper@37950: val gcd=mv_gcd x y; neuper@37950: in neuper@37950: if gcd=[(1,mv_null2(#2(hd(x))))] then neuper@37950: ( neuper@37950: if exists_gcd(x,xs) then new_factors (y::xs @ [x]) neuper@37950: else x::new_factors(y::xs) neuper@37950: ) neuper@37950: else gcd::new_factors(kill_equal(list_div(x::y::xs,gcd))) neuper@37950: end neuper@37950: ) neuper@37950: else neuper@37950: if len=1 then [hd(l)] neuper@37950: else [] neuper@37950: end; neuper@37950: neuper@37950: (*. gets the factors of a list .*) neuper@37950: fun get_factors x = new_factors x; neuper@37950: neuper@37950: (*. multiplies the elements of the list .*) neuper@37950: fun multi_list [] = [] neuper@37950: | multi_list (x::xs) = if xs=[] then x neuper@37950: else mv_mul(x,multi_list xs,LEX_); neuper@37950: neuper@37950: (*. makes a term out of the elements of the list (polynomial representation) .*) neuper@37950: fun make_term ([],vars) = Free(str_of_int 0,HOLogic.realT) neuper@37950: | make_term ((x::xs),vars) = if length(xs)=0 then poly2term(sort (mv_geq LEX_) (x),vars) neuper@37950: else neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2term(sort (mv_geq LEX_) (x),vars) $ neuper@37950: make_term(xs,vars) neuper@37950: ); neuper@37950: neuper@37950: (*. factorizes the denominator (polynomial representation) .*) neuper@37950: fun factorize_den (l,den,vars) = neuper@37950: let neuper@37950: val factor_list=kill_equal( (get_factors l)); neuper@37950: val mlist=multi_list(factor_list); neuper@37950: val (last,rest)=mv_division(den,multi_list(factor_list),LEX_); neuper@37950: in neuper@37950: if rest=[] then neuper@37950: ( neuper@37950: if last=[(1,mv_null2(vars))] then make_term(factor_list,vars) neuper@37950: else make_term(last::factor_list,vars) neuper@37950: ) neuper@37950: else raise error ("RATIONALS_FACTORIZE_DEN_EXCEPTION: Invalid factor by division") neuper@37950: end; neuper@37950: neuper@37950: (*. makes a term out of the elements of the list (expanded polynomial representation) .*) neuper@37950: fun make_exp ([],vars) = Free(str_of_int 0,HOLogic.realT) neuper@37950: | make_exp ((x::xs),vars) = if length(xs)=0 then poly2expanded(sort (mv_geq LEX_) (x),vars) neuper@37950: else neuper@37950: ( neuper@37950: Const ("op *",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: poly2expanded(sort (mv_geq LEX_) (x),vars) $ neuper@37950: make_exp(xs,vars) neuper@37950: ); neuper@37950: neuper@37950: (*. factorizes the denominator (expanded polynomial representation) .*) neuper@37950: fun factorize_den_exp (l,den,vars) = neuper@37950: let neuper@37950: val factor_list=kill_equal( (get_factors l)); neuper@37950: val mlist=multi_list(factor_list); neuper@37950: val (last,rest)=mv_division(den,multi_list(factor_list),LEX_); neuper@37950: in neuper@37950: if rest=[] then neuper@37950: ( neuper@37950: if last=[(1,mv_null2(vars))] then make_exp(factor_list,vars) neuper@37950: else make_exp(last::factor_list,vars) neuper@37950: ) neuper@37950: else raise error ("RATIONALS_FACTORIZE_DEN_EXP_EXCEPTION: Invalid factor by division") neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the common denominator of all elements of the list and multiplies .*) neuper@37950: (*. the nominators and denominators with the correct factor .*) neuper@37950: (*. (polynomial representation) .*) neuper@37950: fun step_add_list_of_fractions []=(Free("0",HOLogic.realT),[]:term list) neuper@37950: | step_add_list_of_fractions [x]= raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXCEPTION: Nothing to add") neuper@37950: | step_add_list_of_fractions (xs) = neuper@37950: let neuper@37950: val den_list=termlist2denominators (xs); (* list of denominators *) neuper@37950: val (denom,var)=calc_lcm(den_list); (* common denominator *) neuper@37950: val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *) neuper@37950: in neuper@37950: com_den(xs,denom,den,var) neuper@37950: end; neuper@37950: neuper@37950: (*. calculates the common denominator of all elements of the list and multiplies .*) neuper@37950: (*. the nominators and denominators with the correct factor .*) neuper@37950: (*. (expanded polynomial representation) .*) neuper@37950: fun step_add_list_of_fractions_exp [] = (Free("0",HOLogic.realT),[]:term list) neuper@37950: | step_add_list_of_fractions_exp [x] = raise error ("RATIONALS_STEP_ADD_LIST_OF_FRACTIONS_EXP_EXCEPTION: Nothing to add") neuper@37950: | step_add_list_of_fractions_exp (xs)= neuper@37950: let neuper@37950: val den_list=termlist2denominators_exp (xs); (* list of denominators *) neuper@37950: val (denom,var)=calc_lcm(den_list); (* common denominator *) neuper@37950: val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *) neuper@37950: in neuper@37950: com_den_exp(xs,denom,den,var) neuper@37950: end; neuper@37950: neuper@37950: (* wird aktuell nicht mehr gebraucht, bei rückänderung schon neuper@37950: ------------------------------------------------------------- neuper@37950: (* WN0210???SK brauch ma des überhaupt *) neuper@37950: fun step_add_list_of_fractions2 []=(Free("0",HOLogic.realT),[]:term list) neuper@37950: | step_add_list_of_fractions2 [x]=(x,[]) neuper@37950: | step_add_list_of_fractions2 (xs) = neuper@37950: let neuper@37950: val den_list=termlist2denominators (xs); (* list of denominators *) neuper@37950: val (denom,var)=calc_lcm(den_list); (* common denominator *) neuper@37950: val den=factorize_den(#1(den_list),denom,var); (* faktorisierter Nenner !!! *) neuper@37950: in neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: com_den2(xs,denom, poly2term(denom,var)(*den*),var) $ neuper@37950: poly2term(denom,var) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: end; neuper@37950: neuper@37950: (* WN0210???SK brauch ma des überhaupt *) neuper@37950: fun step_add_list_of_fractions2_exp []=(Free("0",HOLogic.realT),[]:term list) neuper@37950: | step_add_list_of_fractions2_exp [x]=(x,[]) neuper@37950: | step_add_list_of_fractions2_exp (xs) = neuper@37950: let neuper@37950: val den_list=termlist2denominators_exp (xs); (* list of denominators *) neuper@37950: val (denom,var)=calc_lcm(den_list); (* common denominator *) neuper@37950: val den=factorize_den_exp(#1(den_list),denom,var); (* faktorisierter Nenner !!! *) neuper@37950: in neuper@37950: ( neuper@38014: Const ("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: com_den_exp2(xs,denom, poly2term(denom,var)(*den*),var) $ neuper@37950: poly2expanded(denom,var) neuper@37950: , neuper@37950: [] neuper@37950: ) neuper@37950: end; neuper@37950: ---------------------------------------------- *) neuper@37950: neuper@37950: neuper@37950: (*. converts a term, which contains severel terms seperated by +, into a list of these terms .*) neuper@38014: fun term2list (t as (Const("Rings.inverse_class.divide",_) $ _ $ _)) = [t] neuper@37950: | term2list (t as (Const("Atools.pow",_) $ _ $ _)) = neuper@38014: [Const("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: t $ Free("1",HOLogic.realT) neuper@37950: ] neuper@37950: | term2list (t as (Free(_,_))) = neuper@38014: [Const("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: t $ Free("1",HOLogic.realT) neuper@37950: ] neuper@37950: | term2list (t as (Const("op *",_) $ _ $ _)) = neuper@38014: [Const("Rings.inverse_class.divide",[HOLogic.realT,HOLogic.realT]--->HOLogic.realT) $ neuper@37950: t $ Free("1",HOLogic.realT) neuper@37950: ] neuper@38014: | term2list (Const("Groups.plus_class.plus",_) $ t1 $ t2) = term2list(t1) @ term2list(t2) neuper@38014: | term2list (Const("Groups.minus_class.minus",_) $ t1 $ t2) = neuper@37950: raise error ("RATIONALS_TERM2LIST_EXCEPTION: - not implemented yet") neuper@37950: | term2list _ = raise error ("RATIONALS_TERM2LIST_EXCEPTION: invalid term"); neuper@37950: neuper@37950: (*.factors out the gcd of nominator and denominator: neuper@37950: a/b = (a' * gcd)/(b' * gcd), a,b,gcd are poly[2].*) neuper@37950: fun factout_p_ (thy:theory) t = SOME (step_cancel t,[]:term list); neuper@37950: fun factout_ (thy:theory) t = SOME (step_cancel_expanded t,[]:term list); neuper@37950: neuper@37950: (*.cancels a single fraction with normalform [2] neuper@37950: resulting in a canceled fraction [2], see factout_ .*) neuper@37950: fun cancel_p_ (thy:theory) t = (*WN.2.6.03 no rewrite -> NONE !*) neuper@37950: (let val (t',asm) = direct_cancel(*_expanded ... corrected MG.21.8.03*) t neuper@37950: in if t = t' then NONE else SOME (t',asm) neuper@37950: end) handle _ => NONE; neuper@37950: (*.the same as above with normalform [3] neuper@37950: val cancel_ : neuper@37950: theory -> (*10.02 unused *) neuper@37950: term -> (*fraction in normalform [3] *) neuper@37950: (term * (*fraction in normalform [3] *) neuper@37950: term list) (*casual asumptions in normalform [3] *) neuper@37950: option (*NONE: the function is not applicable *).*) neuper@37950: fun cancel_ (thy:theory) t = SOME (direct_cancel_expanded t) handle _ => NONE; neuper@37950: neuper@37950: (*.transforms sums of at least 2 fractions [3] to neuper@37950: sums with the least common multiple as nominator.*) neuper@37950: fun common_nominator_p_ (thy:theory) t = neuper@37950: ((*writeln("### common_nominator_p_ called");*) neuper@37950: SOME (step_add_list_of_fractions(term2list(t))) handle _ => NONE neuper@37950: ); neuper@37950: fun common_nominator_ (thy:theory) t = neuper@37950: SOME (step_add_list_of_fractions_exp(term2list(t))) handle _ => NONE; neuper@37950: neuper@37950: (*.add 2 or more fractions neuper@37950: val add_fraction_p_ : neuper@37950: theory -> (*10.02 unused *) neuper@37950: term -> (*2 or more fractions with normalform [2] *) neuper@37950: (term * (*one fraction with normalform [2] *) neuper@37950: term list) (*casual assumptions in normalform [2] WN0210???SK *) neuper@37950: option (*NONE: the function is not applicable *).*) neuper@37950: fun add_fraction_p_ (thy:theory) t = neuper@37950: (writeln("### add_fraction_p_ called"); neuper@37950: (let val ts = term2list t neuper@37950: in if 1 < length ts neuper@37950: then SOME (add_list_of_fractions ts) neuper@37950: else NONE (*raise error ("RATIONALS_ADD_EXCEPTION: nothing to add")*) neuper@37950: end) handle _ => NONE neuper@37950: ); neuper@37950: (*.same as add_fraction_p_ but with normalform [3].*) neuper@37950: (*SOME (step_add_list_of_fractions2(term2list(t))); *) neuper@37950: fun add_fraction_ (thy:theory) t = neuper@37950: if length(term2list(t))>1 neuper@37950: then SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE neuper@37950: else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*) neuper@37950: NONE; neuper@37950: fun add_fraction_ (thy:theory) t = neuper@37950: (if 1 < length (term2list t) neuper@37950: then SOME (add_list_of_fractions_exp (term2list t)) neuper@37950: else (*raise error ("RATIONALS_ADD_FRACTION_EXCEPTION: nothing to add")*) neuper@37950: NONE) handle _ => NONE; neuper@37950: neuper@37950: (*SOME (step_add_list_of_fractions2_exp(term2list(t))); *) neuper@37950: neuper@37950: (*. brings the term into a normal form .*) neuper@37950: fun norm_rational_ (thy:theory) t = neuper@37950: SOME (add_list_of_fractions(term2list(t))) handle _ => NONE; neuper@37950: fun norm_expanded_rat_ (thy:theory) t = neuper@37950: SOME (add_list_of_fractions_exp(term2list(t))) handle _ => NONE; neuper@37950: neuper@37950: neuper@37950: (*.evaluates conditions in calculate_Rational.*) neuper@37950: (*make local with FIXX@ME result:term *term list*) neuper@37950: val calc_rat_erls = prep_rls( neuper@37950: Rls {id = "calc_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], (*asm_thm = [], *) neuper@37950: rules = neuper@37950: [Calc ("op =",eval_equal "#equal_"), neuper@37950: Calc ("Atools.is'_const",eval_const "#is_const_"), neuper@37978: Thm ("not_true",num_str @{thm not_true}), neuper@37978: Thm ("not_false",num_str @{thm not_false}) neuper@37950: ], neuper@37950: scr = EmptyScr}); neuper@37950: neuper@37950: neuper@37950: (*.simplifies expressions with numerals; neuper@37950: does NOT rearrange the term by AC-rewriting; thus terms with variables neuper@37950: need to have constants to be commuted together respectively.*) neuper@37950: val calculate_Rational = prep_rls( neuper@37950: merge_rls "calculate_Rational" neuper@37950: (Rls {id = "divide", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = calc_rat_erls, srls = Erls, (*asm_thm = [],*) neuper@37950: calc = [], neuper@37950: rules = neuper@38014: [Calc ("Rings.inverse_class.divide" ,eval_cancel "#divide_e"), neuper@37950: neuper@37965: Thm ("minus_divide_left", neuper@37978: num_str (@{thm minus_divide_left} RS @{thm sym})), neuper@37950: (*SYM - ?x / ?y = - (?x / ?y) may come from subst*) neuper@37950: neuper@37969: Thm ("rat_add",num_str @{thm rat_add}), neuper@37950: (*"[| a is_const; b is_const; c is_const; d is_const |] ==> \ neuper@37950: \"a / c + b / d = (a * d) / (c * d) + (b * c ) / (d * c)"*) neuper@37969: Thm ("rat_add1",num_str @{thm rat_add1}), neuper@37950: (*"[| a is_const; b is_const; c is_const |] ==> \ neuper@37950: \"a / c + b / c = (a + b) / c"*) neuper@37969: Thm ("rat_add2",num_str @{thm rat_add2}), neuper@37950: (*"[| ?a is_const; ?b is_const; ?c is_const |] ==> \ neuper@37950: \?a / ?c + ?b = (?a + ?b * ?c) / ?c"*) neuper@37969: Thm ("rat_add3",num_str @{thm rat_add3}), neuper@37950: (*"[| a is_const; b is_const; c is_const |] ==> \ neuper@37950: \"a + b / c = (a * c) / c + b / c"\ neuper@37950: \.... is_const to be omitted here FIXME*) neuper@37950: neuper@37969: Thm ("rat_mult",num_str @{thm rat_mult}), neuper@37950: (*a / b * (c / d) = a * c / (b * d)*) neuper@37965: Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}), neuper@37950: (*?x * (?y / ?z) = ?x * ?y / ?z*) neuper@37965: Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}), neuper@37950: (*?y / ?z * ?x = ?y * ?x / ?z*) neuper@37950: neuper@37969: Thm ("real_divide_divide1",num_str @{thm real_divide_divide1}), neuper@37950: (*"?y ~= 0 ==> ?u / ?v / (?y / ?z) = ?u / ?v * (?z / ?y)"*) neuper@37965: Thm ("divide_divide_eq_left",num_str @{thm divide_divide_eq_left}), neuper@37950: (*"?x / ?y / ?z = ?x / (?y * ?z)"*) neuper@37950: neuper@37969: Thm ("rat_power", num_str @{thm rat_power}), neuper@37950: (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*) neuper@37950: neuper@37969: Thm ("mult_cross",num_str @{thm mult_cross}), neuper@37950: (*"[| b ~= 0; d ~= 0 |] ==> (a / b = c / d) = (a * d = b * c)*) neuper@37969: Thm ("mult_cross1",num_str @{thm mult_cross1}), neuper@37950: (*" b ~= 0 ==> (a / b = c ) = (a = b * c)*) neuper@37969: Thm ("mult_cross2",num_str @{thm mult_cross2}) neuper@37950: (*" d ~= 0 ==> (a = c / d) = (a * d = c)*) neuper@37950: ], scr = EmptyScr}) neuper@37950: calculate_Poly); neuper@37950: neuper@37950: (*("is_expanded", ("Rational.is'_expanded", eval_is_expanded ""))*) neuper@37950: fun eval_is_expanded (thmid:string) _ neuper@37950: (t as (Const("Rational.is'_expanded", _) $ arg)) thy = neuper@37950: if is_expanded arg neuper@37950: then SOME (mk_thmid thmid "" neuper@37950: ((Syntax.string_of_term (thy2ctxt thy)) arg) "", neuper@37950: Trueprop $ (mk_equality (t, HOLogic.true_const))) neuper@37950: else SOME (mk_thmid thmid "" neuper@37950: ((Syntax.string_of_term (thy2ctxt thy)) arg) "", neuper@37950: Trueprop $ (mk_equality (t, HOLogic.false_const))) neuper@37950: | eval_is_expanded _ _ _ _ = NONE; neuper@37950: neuper@37950: val rational_erls = neuper@37950: merge_rls "rational_erls" calculate_Rational neuper@37950: (append_rls "is_expanded" Atools_erls neuper@37950: [Calc ("Rational.is'_expanded", eval_is_expanded "") neuper@37950: ]); neuper@37950: neuper@37950: neuper@37950: (*.3 'reverse-rewrite-sets' for symbolic computation on rationals: neuper@37950: ================================================================= neuper@37950: A[2] 'cancel_p': . neuper@37950: A[3] 'cancel': . neuper@37950: B[2] 'common_nominator_p': transforms summands in a term [2] neuper@37950: to fractions with the (least) common multiple as nominator. neuper@37950: B[3] 'norm_rational': normalizes arbitrary algebraic terms (without neuper@37950: radicals and transzendental functions) to one canceled fraction, neuper@37950: nominator and denominator in polynomial form. neuper@37950: neuper@37950: In order to meet isac's requirements for interactive and stepwise calculation, neuper@37950: each 'reverse-rewerite-set' consists of an initialization for the interpreter neuper@37950: state and of 4 functions, each of which employs rewriting as much as possible. neuper@37950: The signature of these functions are the same in each 'reverse-rewrite-set' neuper@37950: respectively.*) neuper@37950: neuper@37950: (* ************************************************************************* *) neuper@37950: neuper@37950: local(*. cancel_p neuper@37950: ------------------------ neuper@37950: cancels a single fraction consisting of two (uni- or multivariate) neuper@37950: polynomials WN0609???SK[2] into another such a fraction; examples: neuper@37950: neuper@37950: a^2 + -1*b^2 a + b neuper@37950: -------------------- = --------- neuper@37950: a^2 + -2*a*b + b^2 a + -1*b neuper@37950: neuper@37950: a^2 a neuper@37950: --- = --- neuper@37950: a 1 neuper@37950: neuper@37950: Remark: the reverse ruleset does _NOT_ work properly with other input !.*) neuper@37950: (*WN020824 wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*) neuper@37950: neuper@37950: val {rules, rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "make_polynomial"); neuper@37950: (*WN060829 ... make_deriv does not terminate with 1st expl above, neuper@37950: see rational.sml --- investigate rulesets for cancel_p ---*) neuper@37950: val {rules, rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "rev_rew_p"); neuper@37950: neuper@37950: (*.init_state = fn : term -> istate neuper@37950: initialzies the state of the script interpreter. The state is: neuper@37950: neuper@37950: type rrlsstate = (*state for reverse rewriting*) neuper@37950: (term * (*the current formula*) neuper@37950: term * (*the final term*) neuper@37950: rule list (*'reverse rule list' (#)*) neuper@37950: list * (*may be serveral, eg. in norm_rational*) neuper@37950: (rule * (*Thm (+ Thm generated from Calc) resulting in ...*) neuper@37950: (term * (*... rewrite with ...*) neuper@37950: term list)) (*... assumptions*) neuper@37950: list); (*derivation from given term to normalform neuper@37950: in reverse order with sym_thm; neuper@37950: (#) could be extracted from here by (map #1)*).*) neuper@37950: (* val {rules, rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "rev_rew_p") (*USE ALWAYS, SEE val cancel_p*); neuper@37972: val (thy, eval_rls, ro) =(thy, Atools_erls, ro) (*..val cancel_p*); neuper@37950: val t = t; neuper@37950: *) neuper@37950: fun init_state thy eval_rls ro t = neuper@37950: let val SOME (t',_) = factout_p_ thy t neuper@37950: val SOME (t'',asm) = cancel_p_ thy t neuper@37950: val der = reverse_deriv thy eval_rls rules ro NONE t' neuper@37950: val der = der @ [(Thm ("real_mult_div_cancel2", neuper@37969: num_str @{thm real_mult_div_cancel2}), neuper@37950: (t'',asm))] neuper@37950: val rs = (distinct_Thm o (map #1)) der neuper@37950: val rs = filter_out (eq_Thms ["sym_real_add_zero_left", neuper@37950: "sym_real_mult_0", neuper@37950: "sym_real_mult_1" neuper@37950: (*..insufficient,eg.make_Polynomial*)])rs neuper@37950: in (t,t'',[rs(*here only _ONE_ to ease locate_rule*)],der) end; neuper@37950: neuper@37950: (*.locate_rule = fn : rule list -> term -> rule neuper@37950: -> (rule * (term * term list) option) list. neuper@37950: checks a rule R for being a cancel-rule, and if it is, neuper@37950: then return the list of rules (+ the terms they are rewriting to) neuper@37950: which need to be applied before R should be applied. neuper@37950: precondition: the rule is applicable to the argument-term. neuper@37950: arguments: neuper@37950: rule list: the reverse rule list neuper@37950: -> term : ... to which the rule shall be applied neuper@37950: -> rule : ... to be applied to term neuper@37950: value: neuper@37950: -> (rule : a rule rewriting to ... neuper@37950: * (term : ... the resulting term ... neuper@37950: * term list): ... with the assumptions ( //#0). neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the rule has nothing to do neuper@37950: with cancelation.*) neuper@37950: (* val () = (); neuper@37950: *) neuper@37950: fun locate_rule thy eval_rls ro [rs] t r = neuper@37950: if (id_of_thm r) mem (map (id_of_thm)) rs neuper@37950: then let val ropt = neuper@37950: rewrite_ thy ro eval_rls true (thm_of_thm r) t; neuper@37950: in case ropt of neuper@37950: SOME ta => [(r, ta)] neuper@37950: | NONE => (writeln("### locate_rule: rewrite "^ neuper@37950: (id_of_thm r)^" "^(term2str t)^" = NONE"); neuper@37950: []) end neuper@37950: else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[]) neuper@37950: | locate_rule _ _ _ _ _ _ = neuper@37950: raise error ("locate_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.next_rule = fn : rule list -> term -> rule option neuper@37950: for a given term return the next rules to be done for cancelling. neuper@37950: arguments: neuper@37950: rule list : the reverse rule list neuper@37950: term : the term for which ... neuper@37950: value: neuper@37950: -> rule option: ... this rule is appropriate for cancellation; neuper@37950: there may be no such rule (if the term is canceled already.*) neuper@37972: (* val thy = thy; neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel; neuper@37950: val ([rs],t) = (rss,f); neuper@37950: next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*) neuper@37950: neuper@37972: val (thy, [rs]) = (thy, revsets); neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel; neuper@37950: nex [rs] t; neuper@37950: *) neuper@37950: fun next_rule thy eval_rls ro [rs] t = neuper@37950: let val der = make_deriv thy eval_rls rs ro NONE t; neuper@37950: in case der of neuper@37950: (* val (_,r,_)::_ = der; neuper@37950: *) neuper@37950: (_,r,_)::_ => SOME r neuper@37950: | _ => NONE neuper@37950: end neuper@37950: | next_rule _ _ _ _ _ = neuper@37950: raise error ("next_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.val attach_form = f : rule list -> term -> term neuper@37950: -> (rule * (term * term list)) list neuper@37950: checks an input term TI, if it may belong to a current cancellation, by neuper@37950: trying to derive it from the given term TG. neuper@37950: arguments: neuper@37950: term : TG, the last one in the cancellation agreed upon by user + math-eng neuper@37950: -> term: TI, the next one input by the user neuper@37950: value: neuper@37950: -> (rule : the rule to be applied in order to reach TI neuper@37950: * (term : ... obtained by applying the rule ... neuper@37950: * term list): ... and the respective assumptions. neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the users term does not belong neuper@37950: to a cancellation of the term last agreed upon.*) neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*) neuper@37950: []:(rule * (term * term list)) list; neuper@37950: neuper@37950: in neuper@37950: neuper@37950: val cancel_p = neuper@37950: Rrls {id = "cancel_p", prepat=[], neuper@37950: rew_ord=("ord_make_polynomial", neuper@37972: ord_make_polynomial false thy), neuper@37950: erls = rational_erls, neuper@38014: calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")), neuper@37950: ("TIMES" ,("op *" ,eval_binop "#mult_")), neuper@38014: ("DIVIDE" ,("Rings.inverse_class.divide" ,eval_cancel "#divide_e")), neuper@37950: ("POWER" ,("Atools.pow" ,eval_binop "#power_"))], neuper@37950: (*asm_thm=[("real_mult_div_cancel2","")],*) neuper@37950: scr=Rfuns {init_state = init_state thy Atools_erls ro, neuper@37950: normal_form = cancel_p_ thy, neuper@37950: locate_rule = locate_rule thy Atools_erls ro, neuper@37950: next_rule = next_rule thy Atools_erls ro, neuper@37950: attach_form = attach_form}} neuper@37950: end;(*local*) neuper@37950: neuper@37950: local(*.ad (1) 'cancel' neuper@37950: ------------------------------ neuper@37950: cancels a single fraction consisting of two (uni- or multivariate) neuper@37950: polynomials WN0609???SK[3] into another such a fraction; examples: neuper@37950: neuper@37950: a^2 - b^2 a + b neuper@37950: -------------------- = --------- neuper@37950: a^2 - 2*a*b + b^2 a - *b neuper@37950: neuper@37950: Remark: the reverse ruleset does _NOT_ work properly with other input !.*) neuper@37950: (*WN 24.8.02: wir werden "uberlegen, wie wir ungeeignete inputs zur"uckweisen*) neuper@37950: neuper@37950: (* neuper@37950: val SOME (Rls {rules=rules,rew_ord=(_,ro),...}) = neuper@37950: assoc'(!ruleset',"expand_binoms"); neuper@37950: *) neuper@37950: val {rules=rules,rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "expand_binoms"); neuper@37972: val thy = thy; neuper@37950: neuper@37950: fun init_state thy eval_rls ro t = neuper@37950: let val SOME (t',_) = factout_ thy t; neuper@37950: val SOME (t'',asm) = cancel_ thy t; neuper@37950: val der = reverse_deriv thy eval_rls rules ro NONE t'; neuper@37950: val der = der @ [(Thm ("real_mult_div_cancel2", neuper@37969: num_str @{thm real_mult_div_cancel2}), neuper@37950: (t'',asm))] neuper@37950: val rs = map #1 der; neuper@37950: in (t,t'',[rs],der) end; neuper@37950: neuper@37950: fun locate_rule thy eval_rls ro [rs] t r = neuper@37950: if (id_of_thm r) mem (map (id_of_thm)) rs neuper@37950: then let val ropt = neuper@37950: rewrite_ thy ro eval_rls true (thm_of_thm r) t; neuper@37950: in case ropt of neuper@37950: SOME ta => [(r, ta)] neuper@37950: | NONE => (writeln("### locate_rule: rewrite "^ neuper@37950: (id_of_thm r)^" "^(term2str t)^" = NONE"); neuper@37950: []) end neuper@37950: else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[]) neuper@37950: | locate_rule _ _ _ _ _ _ = neuper@37950: raise error ("locate_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: fun next_rule thy eval_rls ro [rs] t = neuper@37950: let val der = make_deriv thy eval_rls rs ro NONE t; neuper@37950: in case der of neuper@37950: (* val (_,r,_)::_ = der; neuper@37950: *) neuper@37950: (_,r,_)::_ => SOME r neuper@37950: | _ => NONE neuper@37950: end neuper@37950: | next_rule _ _ _ _ _ = neuper@37950: raise error ("next_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*) neuper@37950: []:(rule * (term * term list)) list; neuper@37950: neuper@37979: val pat = parse_patt thy "?r/?s"; neuper@37979: val pre1 = parse_patt thy "?r is_expanded"; neuper@37979: val pre2 = parse_patt thy "?s is_expanded"; neuper@37950: val prepat = [([pre1, pre2], pat)]; neuper@37950: neuper@37950: in neuper@37950: neuper@37950: neuper@37950: val cancel = neuper@37950: Rrls {id = "cancel", prepat=prepat, neuper@37950: rew_ord=("ord_make_polynomial", neuper@37972: ord_make_polynomial false thy), neuper@37950: erls = rational_erls, neuper@38014: calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")), neuper@37950: ("TIMES" ,("op *" ,eval_binop "#mult_")), neuper@38014: ("DIVIDE" ,("Rings.inverse_class.divide" ,eval_cancel "#divide_e")), neuper@37950: ("POWER" ,("Atools.pow" ,eval_binop "#power_"))], neuper@37950: scr=Rfuns {init_state = init_state thy Atools_erls ro, neuper@37950: normal_form = cancel_ thy, neuper@37950: locate_rule = locate_rule thy Atools_erls ro, neuper@37950: next_rule = next_rule thy Atools_erls ro, neuper@37950: attach_form = attach_form}} neuper@37950: end;(*local*) neuper@37950: neuper@37950: local(*.ad [2] 'common_nominator_p' neuper@37950: --------------------------------- neuper@37950: FIXME Beschreibung .*) neuper@37950: neuper@37950: neuper@37950: val {rules=rules,rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "make_polynomial"); neuper@37950: (*WN060829 ... make_deriv does not terminate with 1st expl above, neuper@37950: see rational.sml --- investigate rulesets for cancel_p ---*) neuper@37950: val {rules, rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "rev_rew_p"); neuper@37972: val thy = thy; neuper@37950: neuper@37950: neuper@37950: (*.common_nominator_p_ = fn : theory -> term -> (term * term list) option neuper@37950: as defined above*) neuper@37950: neuper@37950: (*.init_state = fn : term -> istate neuper@37950: initialzies the state of the interactive interpreter. The state is: neuper@37950: neuper@37950: type rrlsstate = (*state for reverse rewriting*) neuper@37950: (term * (*the current formula*) neuper@37950: term * (*the final term*) neuper@37950: rule list (*'reverse rule list' (#)*) neuper@37950: list * (*may be serveral, eg. in norm_rational*) neuper@37950: (rule * (*Thm (+ Thm generated from Calc) resulting in ...*) neuper@37950: (term * (*... rewrite with ...*) neuper@37950: term list)) (*... assumptions*) neuper@37950: list); (*derivation from given term to normalform neuper@37950: in reverse order with sym_thm; neuper@37950: (#) could be extracted from here by (map #1)*).*) neuper@37950: fun init_state thy eval_rls ro t = neuper@37950: let val SOME (t',_) = common_nominator_p_ thy t; neuper@37950: val SOME (t'',asm) = add_fraction_p_ thy t; neuper@37950: val der = reverse_deriv thy eval_rls rules ro NONE t'; neuper@37950: val der = der @ [(Thm ("real_mult_div_cancel2", neuper@37969: num_str @{thm real_mult_div_cancel2}), neuper@37950: (t'',asm))] neuper@37950: val rs = (distinct_Thm o (map #1)) der; neuper@37950: val rs = filter_out (eq_Thms ["sym_real_add_zero_left", neuper@37950: "sym_real_mult_0", neuper@37950: "sym_real_mult_1"]) rs; neuper@37950: in (t,t'',[rs(*here only _ONE_*)],der) end; neuper@37950: neuper@37950: (* use"knowledge/Rational.ML"; neuper@37950: *) neuper@37950: neuper@37950: (*.locate_rule = fn : rule list -> term -> rule neuper@37950: -> (rule * (term * term list) option) list. neuper@37950: checks a rule R for being a cancel-rule, and if it is, neuper@37950: then return the list of rules (+ the terms they are rewriting to) neuper@37950: which need to be applied before R should be applied. neuper@37950: precondition: the rule is applicable to the argument-term. neuper@37950: arguments: neuper@37950: rule list: the reverse rule list neuper@37950: -> term : ... to which the rule shall be applied neuper@37950: -> rule : ... to be applied to term neuper@37950: value: neuper@37950: -> (rule : a rule rewriting to ... neuper@37950: * (term : ... the resulting term ... neuper@37950: * term list): ... with the assumptions ( //#0). neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the rule has nothing to do neuper@37950: with cancelation.*) neuper@37950: (* val () = (); neuper@37950: *) neuper@37950: fun locate_rule thy eval_rls ro [rs] t r = neuper@37950: if (id_of_thm r) mem (map (id_of_thm)) rs neuper@37950: then let val ropt = neuper@37950: rewrite_ thy ro eval_rls true (thm_of_thm r) t; neuper@37950: in case ropt of neuper@37950: SOME ta => [(r, ta)] neuper@37950: | NONE => (writeln("### locate_rule: rewrite "^ neuper@37950: (id_of_thm r)^" "^(term2str t)^" = NONE"); neuper@37950: []) end neuper@37950: else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[]) neuper@37950: | locate_rule _ _ _ _ _ _ = neuper@37950: raise error ("locate_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.next_rule = fn : rule list -> term -> rule option neuper@37950: for a given term return the next rules to be done for cancelling. neuper@37950: arguments: neuper@37950: rule list : the reverse rule list neuper@37950: term : the term for which ... neuper@37950: value: neuper@37950: -> rule option: ... this rule is appropriate for cancellation; neuper@37950: there may be no such rule (if the term is canceled already.*) neuper@37972: (* val thy = thy; neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel; neuper@37950: val ([rs],t) = (rss,f); neuper@37950: next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*) neuper@37950: neuper@37972: val (thy, [rs]) = (thy, revsets); neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel; neuper@37950: nex [rs] t; neuper@37950: *) neuper@37950: fun next_rule thy eval_rls ro [rs] t = neuper@37950: let val der = make_deriv thy eval_rls rs ro NONE t; neuper@37950: in case der of neuper@37950: (* val (_,r,_)::_ = der; neuper@37950: *) neuper@37950: (_,r,_)::_ => SOME r neuper@37950: | _ => NONE neuper@37950: end neuper@37950: | next_rule _ _ _ _ _ = neuper@37950: raise error ("next_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.val attach_form = f : rule list -> term -> term neuper@37950: -> (rule * (term * term list)) list neuper@37950: checks an input term TI, if it may belong to a current cancellation, by neuper@37950: trying to derive it from the given term TG. neuper@37950: arguments: neuper@37950: term : TG, the last one in the cancellation agreed upon by user + math-eng neuper@37950: -> term: TI, the next one input by the user neuper@37950: value: neuper@37950: -> (rule : the rule to be applied in order to reach TI neuper@37950: * (term : ... obtained by applying the rule ... neuper@37950: * term list): ... and the respective assumptions. neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the users term does not belong neuper@37950: to a cancellation of the term last agreed upon.*) neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*) neuper@37950: []:(rule * (term * term list)) list; neuper@37950: neuper@37979: val pat0 = parse_patt thy "?r/?s+?u/?v"; neuper@37979: val pat1 = parse_patt thy "?r/?s+?u "; neuper@37979: val pat2 = parse_patt thy "?r +?u/?v"; neuper@37950: val prepat = [([HOLogic.true_const], pat0), neuper@37950: ([HOLogic.true_const], pat1), neuper@37950: ([HOLogic.true_const], pat2)]; neuper@37950: neuper@37950: in neuper@37950: neuper@37950: (*11.02 schnelle L"osung f"ur RL: Bruch auch gek"urzt; neuper@37950: besser w"are: auf 1 gemeinsamen Bruchstrich, Nenner und Z"ahler unvereinfacht neuper@37950: dh. wie common_nominator_p_, aber auf 1 Bruchstrich*) neuper@37950: val common_nominator_p = neuper@37950: Rrls {id = "common_nominator_p", prepat=prepat, neuper@37950: rew_ord=("ord_make_polynomial", neuper@37972: ord_make_polynomial false thy), neuper@37950: erls = rational_erls, neuper@38014: calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")), neuper@37950: ("TIMES" ,("op *" ,eval_binop "#mult_")), neuper@38014: ("DIVIDE" ,("Rings.inverse_class.divide" ,eval_cancel "#divide_e")), neuper@37950: ("POWER" ,("Atools.pow" ,eval_binop "#power_"))], neuper@37950: scr=Rfuns {init_state = init_state thy Atools_erls ro, neuper@37950: normal_form = add_fraction_p_ thy,(*FIXME.WN0211*) neuper@37950: locate_rule = locate_rule thy Atools_erls ro, neuper@37950: next_rule = next_rule thy Atools_erls ro, neuper@37950: attach_form = attach_form}} neuper@37950: end;(*local*) neuper@37950: neuper@37950: local(*.ad [2] 'common_nominator' neuper@37950: --------------------------------- neuper@37950: FIXME Beschreibung .*) neuper@37950: neuper@37950: neuper@37950: val {rules=rules,rew_ord=(_,ro),...} = neuper@37950: rep_rls (assoc_rls "make_polynomial"); neuper@37972: val thy = thy; neuper@37950: neuper@37950: neuper@37950: (*.common_nominator_ = fn : theory -> term -> (term * term list) option neuper@37950: as defined above*) neuper@37950: neuper@37950: (*.init_state = fn : term -> istate neuper@37950: initialzies the state of the interactive interpreter. The state is: neuper@37950: type rrlsstate = (*state for reverse rewriting*) neuper@37950: (term * (*the current formula*) neuper@37950: term * (*the final term*) neuper@37950: rule list (*'reverse rule list' (#)*) neuper@37950: list * (*may be serveral, eg. in norm_rational*) neuper@37950: (rule * (*Thm (+ Thm generated from Calc) resulting in ...*) neuper@37950: (term * (*... rewrite with ...*) neuper@37950: term list)) (*... assumptions*) neuper@37950: list); (*derivation from given term to normalform neuper@37950: in reverse order with sym_thm; neuper@37950: (#) could be extracted from here by (map #1)*).*) neuper@37950: fun init_state thy eval_rls ro t = neuper@37950: let val SOME (t',_) = common_nominator_ thy t; neuper@37950: val SOME (t'',asm) = add_fraction_ thy t; neuper@37950: val der = reverse_deriv thy eval_rls rules ro NONE t'; neuper@37950: val der = der @ [(Thm ("real_mult_div_cancel2", neuper@37969: num_str @{thm real_mult_div_cancel2}), neuper@37950: (t'',asm))] neuper@37950: val rs = (distinct_Thm o (map #1)) der; neuper@37950: val rs = filter_out (eq_Thms ["sym_real_add_zero_left", neuper@37950: "sym_real_mult_0", neuper@37950: "sym_real_mult_1"]) rs; neuper@37950: in (t,t'',[rs(*here only _ONE_*)],der) end; neuper@37950: neuper@37950: (*.locate_rule = fn : rule list -> term -> rule neuper@37950: -> (rule * (term * term list) option) list. neuper@37950: checks a rule R for being a cancel-rule, and if it is, neuper@37950: then return the list of rules (+ the terms they are rewriting to) neuper@37950: which need to be applied before R should be applied. neuper@37950: precondition: the rule is applicable to the argument-term. neuper@37950: arguments: neuper@37950: rule list: the reverse rule list neuper@37950: -> term : ... to which the rule shall be applied neuper@37950: -> rule : ... to be applied to term neuper@37950: value: neuper@37950: -> (rule : a rule rewriting to ... neuper@37950: * (term : ... the resulting term ... neuper@37950: * term list): ... with the assumptions ( //#0). neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the rule has nothing to do neuper@37950: with cancelation.*) neuper@37950: (* val () = (); neuper@37950: *) neuper@37950: fun locate_rule thy eval_rls ro [rs] t r = neuper@37950: if (id_of_thm r) mem (map (id_of_thm)) rs neuper@37950: then let val ropt = neuper@37950: rewrite_ thy ro eval_rls true (thm_of_thm r) t; neuper@37950: in case ropt of neuper@37950: SOME ta => [(r, ta)] neuper@37950: | NONE => (writeln("### locate_rule: rewrite "^ neuper@37950: (id_of_thm r)^" "^(term2str t)^" = NONE"); neuper@37950: []) end neuper@37950: else (writeln("### locate_rule: "^(id_of_thm r)^" not mem rrls");[]) neuper@37950: | locate_rule _ _ _ _ _ _ = neuper@37950: raise error ("locate_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.next_rule = fn : rule list -> term -> rule option neuper@37950: for a given term return the next rules to be done for cancelling. neuper@37950: arguments: neuper@37950: rule list : the reverse rule list neuper@37950: term : the term for which ... neuper@37950: value: neuper@37950: -> rule option: ... this rule is appropriate for cancellation; neuper@37950: there may be no such rule (if the term is canceled already.*) neuper@37972: (* val thy = thy; neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel; neuper@37950: val ([rs],t) = (rss,f); neuper@37950: next_rule thy eval_rls ro [rs] t;(*eval fun next_rule ... before!*) neuper@37950: neuper@37972: val (thy, [rs]) = (thy, revsets); neuper@37950: val Rrls {rew_ord=(_,ro),...} = cancel_p; neuper@37950: nex [rs] t; neuper@37950: *) neuper@37950: fun next_rule thy eval_rls ro [rs] t = neuper@37950: let val der = make_deriv thy eval_rls rs ro NONE t; neuper@37950: in case der of neuper@37950: (* val (_,r,_)::_ = der; neuper@37950: *) neuper@37950: (_,r,_)::_ => SOME r neuper@37950: | _ => NONE neuper@37950: end neuper@37950: | next_rule _ _ _ _ _ = neuper@37950: raise error ("next_rule: doesnt match rev-sets in istate"); neuper@37950: neuper@37950: (*.val attach_form = f : rule list -> term -> term neuper@37950: -> (rule * (term * term list)) list neuper@37950: checks an input term TI, if it may belong to a current cancellation, by neuper@37950: trying to derive it from the given term TG. neuper@37950: arguments: neuper@37950: term : TG, the last one in the cancellation agreed upon by user + math-eng neuper@37950: -> term: TI, the next one input by the user neuper@37950: value: neuper@37950: -> (rule : the rule to be applied in order to reach TI neuper@37950: * (term : ... obtained by applying the rule ... neuper@37950: * term list): ... and the respective assumptions. neuper@37950: ) list : there may be several such rules; neuper@37950: the list is empty, if the users term does not belong neuper@37950: to a cancellation of the term last agreed upon.*) neuper@37950: fun attach_form (_:rule list list) (_:term) (_:term) = (*still missing*) neuper@37950: []:(rule * (term * term list)) list; neuper@37950: neuper@37979: val pat0 = parse_patt thy "?r/?s+?u/?v"; neuper@37979: val pat01 = parse_patt thy "?r/?s-?u/?v"; neuper@37979: val pat1 = parse_patt thy "?r/?s+?u "; neuper@37979: val pat11 = parse_patt thy "?r/?s-?u "; neuper@37979: val pat2 = parse_patt thy "?r +?u/?v"; neuper@37979: val pat21 = parse_patt thy "?r -?u/?v"; neuper@37950: val prepat = [([HOLogic.true_const], pat0), neuper@37950: ([HOLogic.true_const], pat01), neuper@37950: ([HOLogic.true_const], pat1), neuper@37950: ([HOLogic.true_const], pat11), neuper@37950: ([HOLogic.true_const], pat2), neuper@37950: ([HOLogic.true_const], pat21)]; neuper@37950: neuper@37950: neuper@37950: in neuper@37950: neuper@37950: val common_nominator = neuper@37950: Rrls {id = "common_nominator", prepat=prepat, neuper@37950: rew_ord=("ord_make_polynomial", neuper@37972: ord_make_polynomial false thy), neuper@37950: erls = rational_erls, neuper@38014: calc = [("PLUS" ,("Groups.plus_class.plus" ,eval_binop "#add_")), neuper@37950: ("TIMES" ,("op *" ,eval_binop "#mult_")), neuper@38014: ("DIVIDE" ,("Rings.inverse_class.divide" ,eval_cancel "#divide_e")), neuper@37950: ("POWER" ,("Atools.pow" ,eval_binop "#power_"))], neuper@37950: (*asm_thm=[("real_mult_div_cancel2","")],*) neuper@37950: scr=Rfuns {init_state = init_state thy Atools_erls ro, neuper@37950: normal_form = add_fraction_ (*NOT common_nominator_*) thy, neuper@37950: locate_rule = locate_rule thy Atools_erls ro, neuper@37950: next_rule = next_rule thy Atools_erls ro, neuper@37950: attach_form = attach_form}} neuper@37950: neuper@37950: end;(*local*) neuper@37950: end;(*struct*) neuper@37950: neuper@37950: open RationalI; neuper@37950: (*##*) neuper@37950: neuper@37950: (*.the expression contains + - * ^ / only ?.*) neuper@37950: fun is_ratpolyexp (Free _) = true neuper@38014: | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ Free _ $ Free _) = true neuper@38014: | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ Free _ $ Free _) = true neuper@37950: | is_ratpolyexp (Const ("op *",_) $ Free _ $ Free _) = true neuper@37950: | is_ratpolyexp (Const ("Atools.pow",_) $ Free _ $ Free _) = true neuper@38014: | is_ratpolyexp (Const ("Rings.inverse_class.divide",_) $ Free _ $ Free _) = true neuper@38014: | is_ratpolyexp (Const ("Groups.plus_class.plus",_) $ t1 $ t2) = neuper@37950: ((is_ratpolyexp t1) andalso (is_ratpolyexp t2)) neuper@38014: | is_ratpolyexp (Const ("Groups.minus_class.minus",_) $ t1 $ t2) = neuper@37950: ((is_ratpolyexp t1) andalso (is_ratpolyexp t2)) neuper@37950: | is_ratpolyexp (Const ("op *",_) $ t1 $ t2) = neuper@37950: ((is_ratpolyexp t1) andalso (is_ratpolyexp t2)) neuper@37950: | is_ratpolyexp (Const ("Atools.pow",_) $ t1 $ t2) = neuper@37950: ((is_ratpolyexp t1) andalso (is_ratpolyexp t2)) neuper@38014: | is_ratpolyexp (Const ("Rings.inverse_class.divide",_) $ t1 $ t2) = neuper@37950: ((is_ratpolyexp t1) andalso (is_ratpolyexp t2)) neuper@37950: | is_ratpolyexp _ = false; neuper@37950: neuper@37950: (*("is_ratpolyexp", ("Rational.is'_ratpolyexp", eval_is_ratpolyexp ""))*) neuper@37950: fun eval_is_ratpolyexp (thmid:string) _ neuper@37950: (t as (Const("Rational.is'_ratpolyexp", _) $ arg)) thy = neuper@37950: if is_ratpolyexp arg neuper@37950: then SOME (mk_thmid thmid "" neuper@37950: ((Syntax.string_of_term (thy2ctxt thy)) arg) "", neuper@37950: Trueprop $ (mk_equality (t, HOLogic.true_const))) neuper@37950: else SOME (mk_thmid thmid "" neuper@37950: ((Syntax.string_of_term (thy2ctxt thy)) arg) "", neuper@37950: Trueprop $ (mk_equality (t, HOLogic.false_const))) neuper@37950: | eval_is_ratpolyexp _ _ _ _ = NONE; neuper@37950: neuper@37950: neuper@37950: neuper@37950: (*-------------------18.3.03 --> struct <-----------vvv--*) neuper@37950: val add_fractions_p = common_nominator_p; (*FIXXXME:eilig f"ur norm_Rational*) neuper@37950: neuper@37980: (*WN100906 removed in favour of discard_minus = discard_minus_ formerly neuper@37980: discard binary minus, shift unary minus into -1*; neuper@37950: unary minus before numerals are put into the numeral by parsing; neuper@37980: contains absolute minimum of thms for context in norm_Rational neuper@37950: val discard_minus = prep_rls( neuper@37950: Rls {id = "discard_minus", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37980: erls = e_rls, srls = Erls, calc = [], neuper@37969: rules = [Thm ("real_diff_minus", num_str @{thm real_diff_minus}), neuper@37950: (*"a - b = a + -1 * b"*) neuper@37969: Thm ("sym_real_mult_minus1", neuper@37969: num_str (@{thm real_mult_minus1} RS @{thm sym})) neuper@37950: (*- ?z = "-1 * ?z"*) neuper@37950: ], neuper@37979: scr = EmptyScr neuper@37950: }):rls; neuper@37980: *) neuper@37980: neuper@37950: (*erls for calculate_Rational; make local with FIXX@ME result:term *term list*) neuper@37950: val powers_erls = prep_rls( neuper@37950: Rls {id = "powers_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*) neuper@37950: rules = [Calc ("Atools.is'_atom",eval_is_atom "#is_atom_"), neuper@37950: Calc ("Atools.is'_even",eval_is_even "#is_even_"), neuper@37950: Calc ("op <",eval_equ "#less_"), neuper@37979: Thm ("not_false", num_str @{thm not_false}), neuper@37979: Thm ("not_true", num_str @{thm not_true}), neuper@38014: Calc ("Groups.plus_class.plus",eval_binop "#add_") neuper@37950: ], neuper@37979: scr = EmptyScr neuper@37950: }:rls); neuper@37950: (*.all powers over + distributed; atoms over * collected, other distributed neuper@37950: contains absolute minimum of thms for context in norm_Rational .*) neuper@37950: val powers = prep_rls( neuper@37950: Rls {id = "powers", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = powers_erls, srls = Erls, calc = [], (*asm_thm = [],*) neuper@37969: rules = [Thm ("realpow_multI", num_str @{thm realpow_multI}), neuper@37950: (*"(r * s) ^^^ n = r ^^^ n * s ^^^ n"*) neuper@37969: Thm ("realpow_pow",num_str @{thm realpow_pow}), neuper@37950: (*"(a ^^^ b) ^^^ c = a ^^^ (b * c)"*) neuper@37969: Thm ("realpow_oneI",num_str @{thm realpow_oneI}), neuper@37950: (*"r ^^^ 1 = r"*) neuper@37969: Thm ("realpow_minus_even",num_str @{thm realpow_minus_even}), neuper@37950: (*"n is_even ==> (- r) ^^^ n = r ^^^ n" ?-->discard_minus?*) neuper@37969: Thm ("realpow_minus_odd",num_str @{thm realpow_minus_odd}), neuper@37950: (*"Not (n is_even) ==> (- r) ^^^ n = -1 * r ^^^ n"*) neuper@37950: neuper@37950: (*----- collect atoms over * -----*) neuper@37969: Thm ("realpow_two_atom",num_str @{thm realpow_two_atom}), neuper@37950: (*"r is_atom ==> r * r = r ^^^ 2"*) neuper@37969: Thm ("realpow_plus_1",num_str @{thm realpow_plus_1}), neuper@37950: (*"r is_atom ==> r * r ^^^ n = r ^^^ (n + 1)"*) neuper@37969: Thm ("realpow_addI_atom",num_str @{thm realpow_addI_atom}), neuper@37950: (*"r is_atom ==> r ^^^ n * r ^^^ m = r ^^^ (n + m)"*) neuper@37950: neuper@37950: (*----- distribute none-atoms -----*) neuper@37969: Thm ("realpow_def_atom",num_str @{thm realpow_def_atom}), neuper@37950: (*"[| 1 < n; not(r is_atom) |]==>r ^^^ n = r * r ^^^ (n + -1)"*) neuper@37969: Thm ("realpow_eq_oneI",num_str @{thm realpow_eq_oneI}), neuper@37950: (*"1 ^^^ n = 1"*) neuper@38014: Calc ("Groups.plus_class.plus",eval_binop "#add_") neuper@37950: ], neuper@37979: scr = EmptyScr neuper@37950: }:rls); neuper@37950: (*.contains absolute minimum of thms for context in norm_Rational.*) neuper@37950: val rat_mult_divide = prep_rls( neuper@37950: Rls {id = "rat_mult_divide", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*) neuper@37969: rules = [Thm ("rat_mult",num_str @{thm rat_mult}), neuper@37950: (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*) neuper@37965: Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}), neuper@37950: (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2], neuper@37950: otherwise inv.to a / b / c = ...*) neuper@37965: Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}), neuper@37950: (*"?a / ?b * ?c = ?a * ?c / ?b" order weights x^^^n too much neuper@37950: and does not commute a / b * c ^^^ 2 !*) neuper@37950: neuper@37979: Thm ("divide_divide_eq_right", neuper@37979: num_str @{thm divide_divide_eq_right}), neuper@37950: (*"?x / (?y / ?z) = ?x * ?z / ?y"*) neuper@37979: Thm ("divide_divide_eq_left", neuper@37979: num_str @{thm divide_divide_eq_left}), neuper@37950: (*"?x / ?y / ?z = ?x / (?y * ?z)"*) neuper@38014: Calc ("Rings.inverse_class.divide" ,eval_cancel "#divide_e") neuper@37950: ], neuper@37979: scr = EmptyScr neuper@37950: }:rls); neuper@37979: neuper@37950: (*.contains absolute minimum of thms for context in norm_Rational.*) neuper@37950: val reduce_0_1_2 = prep_rls( neuper@37950: Rls{id = "reduce_0_1_2", preconds = [], rew_ord = ("dummy_ord", dummy_ord), neuper@37950: erls = e_rls,srls = Erls,calc = [],(*asm_thm = [],*) neuper@37965: rules = [(*Thm ("divide_1",num_str @{thm divide_1}), neuper@37950: "?x / 1 = ?x" unnecess.for normalform*) neuper@37965: Thm ("mult_1_left",num_str @{thm mult_1_left}), neuper@37950: (*"1 * z = z"*) neuper@37969: (*Thm ("real_mult_minus1",num_str @{thm real_mult_minus1}), neuper@37950: "-1 * z = - z"*) neuper@37969: (*Thm ("real_minus_mult_cancel",num_str @{thm real_minus_mult_cancel}), neuper@37950: "- ?x * - ?y = ?x * ?y"*) neuper@37950: neuper@37965: Thm ("mult_zero_left",num_str @{thm mult_zero_left}), neuper@37950: (*"0 * z = 0"*) neuper@37965: Thm ("add_0_left",num_str @{thm add_0_left}), neuper@37950: (*"0 + z = z"*) neuper@37965: (*Thm ("right_minus",num_str @{thm right_minus}), neuper@37950: "?z + - ?z = 0"*) neuper@37950: neuper@37969: Thm ("sym_real_mult_2", neuper@37969: num_str (@{thm real_mult_2} RS @{thm sym})), neuper@37950: (*"z1 + z1 = 2 * z1"*) neuper@37969: Thm ("real_mult_2_assoc",num_str @{thm real_mult_2_assoc}), neuper@37950: (*"z1 + (z1 + k) = 2 * z1 + k"*) neuper@37950: neuper@37965: Thm ("divide_zero_left",num_str @{thm divide_zero_left}) neuper@37950: (*"0 / ?x = 0"*) neuper@37950: ], scr = EmptyScr}:rls); neuper@37950: neuper@37950: (*erls for calculate_Rational; neuper@37950: make local with FIXX@ME result:term *term list WN0609???SKMG*) neuper@37950: val norm_rat_erls = prep_rls( neuper@37950: Rls {id = "norm_rat_erls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], (*asm_thm = [],*) neuper@37950: rules = [Calc ("Atools.is'_const",eval_const "#is_const_") neuper@37979: ], scr = EmptyScr}:rls); neuper@37979: neuper@37950: (*.consists of rls containing the absolute minimum of thms.*) neuper@37950: (*040209: this version has been used by RL for his equations, neuper@37950: which is now replaced by MGs version below neuper@37950: vvv OLD VERSION !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!*) neuper@37950: val norm_Rational = prep_rls( neuper@37950: Rls {id = "norm_Rational", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@37979: erls = norm_rat_erls, srls = Erls, calc = [], neuper@37950: rules = [(*sequence given by operator precedence*) neuper@37950: Rls_ discard_minus, neuper@37950: Rls_ powers, neuper@37950: Rls_ rat_mult_divide, neuper@37950: Rls_ expand, neuper@37950: Rls_ reduce_0_1_2, neuper@37950: (*^^^^^^^^^ from RL -- not the latest one vvvvvvvvv*) neuper@37950: Rls_ order_add_mult, neuper@37950: Rls_ collect_numerals, neuper@37950: Rls_ add_fractions_p, neuper@37950: Rls_ cancel_p neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37979: neuper@37950: val norm_Rational_parenthesized = prep_rls( neuper@37950: Seq {id = "norm_Rational_parenthesized", preconds = []:term list, neuper@37950: rew_ord = ("dummy_ord", dummy_ord), neuper@37950: erls = Atools_erls, srls = Erls, neuper@37979: calc = [], neuper@37950: rules = [Rls_ norm_Rational, (*from RL -- not the latest one*) neuper@37950: Rls_ discard_parentheses neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37950: neuper@37950: (*WN030318???SK: simplifies all but cancel and common_nominator*) neuper@37950: val simplify_rational = neuper@37950: merge_rls "simplify_rational" expand_binoms neuper@37950: (append_rls "divide" calculate_Rational neuper@37965: [Thm ("divide_1",num_str @{thm divide_1}), neuper@37950: (*"?x / 1 = ?x"*) neuper@37978: Thm ("rat_mult",num_str @{thm rat_mult}), neuper@37950: (*(1)"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*) neuper@37965: Thm ("times_divide_eq_right",num_str @{thm times_divide_eq_right}), neuper@37950: (*(2)"?a * (?c / ?d) = ?a * ?c / ?d" must be [2], neuper@37950: otherwise inv.to a / b / c = ...*) neuper@37965: Thm ("times_divide_eq_left",num_str @{thm times_divide_eq_left}), neuper@37950: (*"?a / ?b * ?c = ?a * ?c / ?b"*) neuper@37969: Thm ("add_minus",num_str @{thm add_minus}), neuper@37950: (*"?a + ?b - ?b = ?a"*) neuper@37969: Thm ("add_minus1",num_str @{thm add_minus1}), neuper@37950: (*"?a - ?b + ?b = ?a"*) neuper@37978: Thm ("divide_minus1",num_str @{thm divide_minus1}) neuper@37950: (*"?x / -1 = - ?x"*) neuper@37950: (* neuper@37950: , neuper@37969: Thm ("",num_str @{thm }) neuper@37950: *) neuper@37950: ]); neuper@37950: neuper@37950: (*---------vvv-------------MG ab 1.07.2003--------------vvv-----------*) neuper@37950: neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: (* Simplifier für beliebige Buchterme *) neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: (*----------------------- norm_Rational_mg ---------------------------*) neuper@37950: (*. description of the simplifier see MG-DA.p.56ff .*) neuper@37950: (* ------------------------------------------------------------------- *) neuper@37950: val common_nominator_p_rls = prep_rls( neuper@37950: Rls {id = "common_nominator_p_rls", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], neuper@37950: rules = neuper@37950: [Rls_ common_nominator_p neuper@37950: (*FIXME.WN0401 ? redesign Rrls - use exhaustively on a term ? neuper@37950: FIXME.WN0510 unnecessary nesting: introduce RRls_ : rls -> rule*) neuper@37950: ], neuper@37950: scr = EmptyScr}); neuper@37950: (* ------------------------------------------------------------------- *) neuper@37950: val cancel_p_rls = prep_rls( neuper@37950: Rls {id = "cancel_p_rls", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], neuper@37950: rules = neuper@37950: [Rls_ cancel_p neuper@37950: (*FIXME.WN.0401 ? redesign Rrls - use exhaustively on a term ?*) neuper@37950: ], neuper@37950: scr = EmptyScr}); neuper@37950: (* -------------------------------------------------------------------- *) neuper@37950: (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions; neuper@37950: used in initial part norm_Rational_mg, see example DA-M02-main.p.60.*) neuper@37950: val rat_mult_poly = prep_rls( neuper@37950: Rls {id = "rat_mult_poly", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = append_rls "e_rls-is_polyexp" e_rls neuper@37950: [Calc ("Poly.is'_polyexp", eval_is_polyexp "")], neuper@37950: srls = Erls, calc = [], neuper@37950: rules = neuper@37969: [Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}), neuper@37950: (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*) neuper@37969: Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}) neuper@37950: (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) neuper@37950: ], neuper@37950: scr = EmptyScr}); neuper@37979: neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: (*. makes 'normal' fractions; 'is_polyexp' inhibits double fractions; neuper@37950: used in looping part norm_Rational_rls, see example DA-M02-main.p.60 neuper@37950: .. WHERE THE LATTER DOES ALWAYS WORK, BECAUSE erls = e_rls, neuper@37950: I.E. THE RESPECTIVE ASSUMPTION IS STORED AND Thm APPLIED; WN051028 neuper@37950: ... WN0609???MG.*) neuper@37950: val rat_mult_div_pow = prep_rls( neuper@37950: Rls {id = "rat_mult_div_pow", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, neuper@37950: (*FIXME.WN051028 append_rls "e_rls-is_polyexp" e_rls neuper@37950: [Calc ("Poly.is'_polyexp", eval_is_polyexp "")], neuper@37950: with this correction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ we get neuper@37950: error "rational.sml.sml: diff.behav. in norm_Rational_mg 29" etc. neuper@37950: thus we decided to go on with this flaw*) neuper@37950: srls = Erls, calc = [], neuper@37969: rules = [Thm ("rat_mult",num_str @{thm rat_mult}), neuper@37950: (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*) neuper@37969: Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}), neuper@37950: (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*) neuper@37969: Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}), neuper@37950: (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*) neuper@37950: neuper@37979: Thm ("real_divide_divide1_mg", neuper@37979: num_str @{thm real_divide_divide1_mg}), neuper@37950: (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*) neuper@37979: Thm ("divide_divide_eq_right", neuper@37979: num_str @{thm divide_divide_eq_right}), neuper@37950: (*"?x / (?y / ?z) = ?x * ?z / ?y"*) neuper@37979: Thm ("divide_divide_eq_left", neuper@37979: num_str @{thm divide_divide_eq_left}), neuper@37950: (*"?x / ?y / ?z = ?x / (?y * ?z)"*) neuper@38014: Calc ("Rings.inverse_class.divide" ,eval_cancel "#divide_e"), neuper@37950: neuper@37969: Thm ("rat_power", num_str @{thm rat_power}) neuper@37950: (*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*) neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: val rat_reduce_1 = prep_rls( neuper@37950: Rls {id = "rat_reduce_1", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = e_rls, srls = Erls, calc = [], neuper@37965: rules = [Thm ("divide_1",num_str @{thm divide_1}), neuper@37950: (*"?x / 1 = ?x"*) neuper@37965: Thm ("mult_1_left",num_str @{thm mult_1_left}) neuper@37950: (*"1 * z = z"*) neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: (*. looping part of norm_Rational(*_mg*) .*) neuper@37950: val norm_Rational_rls = prep_rls( neuper@37950: Rls {id = "norm_Rational_rls", preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = norm_rat_erls, srls = Erls, calc = [], neuper@37950: rules = [Rls_ common_nominator_p_rls, neuper@37950: Rls_ rat_mult_div_pow, neuper@37950: Rls_ make_rat_poly_with_parentheses, neuper@37950: Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*) neuper@37950: Rls_ rat_reduce_1 neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: (*040109 'norm_Rational'(by RL) replaced by 'norm_Rational_mg'(MG) neuper@37950: just be renaming:*) neuper@37950: val norm_Rational(*_mg*) = prep_rls( neuper@37950: Seq {id = "norm_Rational"(*_mg*), preconds = [], neuper@37950: rew_ord = ("dummy_ord",dummy_ord), neuper@37950: erls = norm_rat_erls, srls = Erls, calc = [], neuper@37980: rules = [Rls_ discard_minus, neuper@37950: Rls_ rat_mult_poly,(* removes double fractions like a/b/c *) neuper@37950: Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*) neuper@37950: Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*) neuper@37950: Rls_ norm_Rational_rls, (* the main rls, looping (#) *) neuper@37979: Rls_ discard_parentheses1 (* mult only *) neuper@37950: ], neuper@37979: scr = EmptyScr}:rls); neuper@37950: (* ------------------------------------------------------------------ *) neuper@37950: neuper@37967: ruleset' := overwritelthy @{theory} (!ruleset', neuper@37950: [("calculate_Rational", calculate_Rational), neuper@37950: ("calc_rat_erls",calc_rat_erls), neuper@37950: ("rational_erls", rational_erls), neuper@37950: ("cancel_p", cancel_p), neuper@37950: ("cancel", cancel), neuper@37950: ("common_nominator_p", common_nominator_p), neuper@37950: ("common_nominator_p_rls", common_nominator_p_rls), neuper@37950: ("common_nominator" , common_nominator), neuper@37950: ("discard_minus", discard_minus), neuper@37950: ("powers_erls", powers_erls), neuper@37950: ("powers", powers), neuper@37950: ("rat_mult_divide", rat_mult_divide), neuper@37950: ("reduce_0_1_2", reduce_0_1_2), neuper@37950: ("rat_reduce_1", rat_reduce_1), neuper@37950: ("norm_rat_erls", norm_rat_erls), neuper@37950: ("norm_Rational", norm_Rational), neuper@37950: ("norm_Rational_rls", norm_Rational_rls), neuper@37950: ("norm_Rational_parenthesized", norm_Rational_parenthesized), neuper@37950: ("rat_mult_poly", rat_mult_poly), neuper@37950: ("rat_mult_div_pow", rat_mult_div_pow), neuper@37950: ("cancel_p_rls", cancel_p_rls) neuper@37950: ]); neuper@37950: neuper@37950: calclist':= overwritel (!calclist', neuper@37950: [("is_expanded", ("Rational.is'_expanded", eval_is_expanded "")) neuper@37950: ]); neuper@37950: neuper@37950: (** problems **) neuper@37950: neuper@37950: store_pbt neuper@37972: (prep_pbt thy "pbl_simp_rat" [] e_pblID neuper@37950: (["rational","simplification"], neuper@37979: [("#Given" ,["TERM t_t"]), neuper@37979: ("#Where" ,["t_t is_ratpolyexp"]), neuper@37979: ("#Find" ,["normalform n_n"]) neuper@37950: ], neuper@37950: append_rls "e_rls" e_rls [(*for preds in where_*)], neuper@37979: SOME "Simplifyt_t", neuper@37950: [["simplification","of_rationals"]])); neuper@37950: neuper@37950: (** methods **) neuper@37979: val --------------------------------------------------+++++-+++++ = "00000"; neuper@37950: neuper@37950: (*WN061025 this methods script is copied from (auto-generated) script neuper@37950: of norm_Rational in order to ease repair on inform*) neuper@37950: store_met neuper@37972: (prep_met thy "met_simp_rat" [] e_metID neuper@37950: (["simplification","of_rationals"], neuper@37979: [("#Given" ,["TERM t_t"]), neuper@37979: ("#Where" ,["t_t is_ratpolyexp"]), neuper@37979: ("#Find" ,["normalform n_n"]) neuper@37950: ], neuper@37950: {rew_ord'="tless_true", neuper@37950: rls' = e_rls, neuper@37950: calc = [], srls = e_rls, neuper@37950: prls = append_rls "simplification_of_rationals_prls" e_rls neuper@37950: [(*for preds in where_*) neuper@37950: Calc ("Rational.is'_ratpolyexp", neuper@37950: eval_is_ratpolyexp "")], neuper@37950: crls = e_rls, nrls = norm_Rational_rls}, neuper@37979: "Script SimplifyScript (t_t::real) = " ^ neuper@37980: " ((Try (Rewrite_Set discard_minus False) @@ " ^ neuper@37950: " Try (Rewrite_Set rat_mult_poly False) @@ " ^ neuper@37950: " Try (Rewrite_Set make_rat_poly_with_parentheses False) @@ " ^ neuper@37950: " Try (Rewrite_Set cancel_p_rls False) @@ " ^ neuper@37950: " (Repeat " ^ neuper@37950: " ((Try (Rewrite_Set common_nominator_p_rls False) @@ " ^ neuper@37950: " Try (Rewrite_Set rat_mult_div_pow False) @@ " ^ neuper@37950: " Try (Rewrite_Set make_rat_poly_with_parentheses False) @@" ^ neuper@37950: " Try (Rewrite_Set cancel_p_rls False) @@ " ^ neuper@37950: " Try (Rewrite_Set rat_reduce_1 False)))) @@ " ^ neuper@37979: " Try (Rewrite_Set discard_parentheses1 False)) " ^ neuper@37979: " t_t)" neuper@37950: )); neuper@37979: neuper@37950: *} neuper@37950: neuper@37950: end