neuper@42376: (* Partial_Fractions jan@42344: author: Jan Rocnik, isac team jan@42344: Copyright (c) isac team 2011 jan@42344: Use is subject to license terms. jan@42344: jan@42344: 1234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 jan@42344: 10 20 30 40 50 60 70 80 90 100 jan@42344: *) neuper@42376: header {* Partial Fraction Decomposition *} jan@42344: neuper@42289: neuper@42376: theory Partial_Fractions imports RootRatEq begin jan@42353: jan@42353: ML{* neuper@42376: (* jan@42353: signature PARTFRAC = jan@42353: sig jan@42353: val ansatz_rls : rls jan@42353: val ansatz_rls_ : theory -> term -> (term * term list) option jan@42353: end neuper@42376: *) jan@42353: *} jan@42353: neuper@42376: subsection {* eval_ functions *} jan@42344: consts neuper@42359: factors_from_solution :: "bool list => real" neuper@42359: drop_questionmarks :: "'a => 'a" jan@42344: neuper@42376: text {* these might be used for variants of fac_from_sol *} neuper@42289: ML {* jan@42344: fun mk_minus_1 T = Free("-1", T); (*TODO DELETE WITH numbers_to_string*) jan@42344: fun flip_sign t = (*TODO improve for use in factors_from_solution: -(-1) etc*) jan@42344: let val minus_1 = t |> type_of |> mk_minus_1 jan@42344: in HOLogic.mk_binop "Groups.times_class.times" (minus_1, t) end; neuper@42376: *} neuper@42376: neuper@42376: text {* from solutions (e.g. [z = 1, z = -2]) make linear factors (e.g. (z - 1)*(z - -2)) *} neuper@42376: ML {* jan@42344: fun fac_from_sol s = jan@42344: let val (lhs, rhs) = HOLogic.dest_eq s jan@42367: in HOLogic.mk_binop "Groups.minus_class.minus" (lhs, rhs) end; jan@42344: jan@42344: fun mk_prod prod [] = jan@42344: if prod = e_term then error "mk_prod called with []" else prod jan@42344: | mk_prod prod (t :: []) = jan@42344: if prod = e_term then t else HOLogic.mk_binop "Groups.times_class.times" (prod, t) jan@42344: | mk_prod prod (t1 :: t2 :: ts) = jan@42344: if prod = e_term jan@42344: then jan@42344: let val p = HOLogic.mk_binop "Groups.times_class.times" (t1, t2) jan@42344: in mk_prod p ts end jan@42344: else jan@42344: let val p = HOLogic.mk_binop "Groups.times_class.times" (prod, t1) jan@42344: in mk_prod p (t2 :: ts) end jan@42344: jan@42344: fun factors_from_solution sol = jan@42344: let val ts = HOLogic.dest_list sol jan@42344: in mk_prod e_term (map fac_from_sol ts) end; jan@42344: neuper@42376: (*("factors_from_solution", ("Partial_Fractions.factors_from_solution", neuper@42376: eval_factors_from_solution ""))*) jan@42352: fun eval_factors_from_solution (thmid:string) _ jan@42352: (t as Const ("Partial_Fractions.factors_from_solution", _) $ sol) thy = jan@42352: ((let val prod = factors_from_solution sol neuper@52070: in SOME (mk_thmid thmid "" (term_to_string''' thy prod) "", neuper@52070: Trueprop $ (mk_equality (t, prod))) jan@42352: end) jan@42352: handle _ => NONE) jan@42352: | eval_factors_from_solution _ _ _ _ = NONE; jan@42344: *} jan@42344: neuper@42376: text {* 'ansatz' introduces '?Vars' (questionable design); drop these again *} neuper@42359: ML {* neuper@42359: (*("drop_questionmarks", ("Partial_Fractions.drop_questionmarks", eval_drop_questionmarks ""))*) neuper@42359: fun eval_drop_questionmarks (thmid:string) _ neuper@42359: (t as Const ("Partial_Fractions.drop_questionmarks", _) $ tm) thy = neuper@42359: if contains_Var tm neuper@42359: then neuper@42359: let neuper@42359: val tm' = var2free tm neuper@52070: in SOME (mk_thmid thmid "" (term_to_string''' thy tm') "", neuper@52070: Trueprop $ (mk_equality (t, tm'))) neuper@42359: end neuper@42359: else NONE neuper@42359: | eval_drop_questionmarks _ _ _ _ = NONE; neuper@42376: *} neuper@42359: neuper@42376: text {* store eval_ functions for calls from Scripts *} s1210629013@52145: setup {* KEStore_Elems.add_calcs s1210629013@52145: [("drop_questionmarks", ("Partial_Fractions.drop'_questionmarks", eval_drop_questionmarks ""))] *} neuper@42359: neuper@42376: subsection {* 'ansatz' for partial fractions *} jan@42353: axiomatization where jan@42358: ansatz_2nd_order: "n / (a*b) = A/a + B/b" and neuper@42376: ansatz_3rd_order: "n / (a*b*c) = A/a + B/b + C/c" and neuper@42376: ansatz_4th_order: "n / (a*b*c*d) = A/a + B/b + C/c + D/d" and neuper@42386: (*version 1*) neuper@42376: equival_trans_2nd_order: "(n/(a*b) = A/a + B/b) = (n = A*b + B*a)" and neuper@42376: equival_trans_3rd_order: "(n/(a*b*c) = A/a + B/b + C/c) = (n = A*b*c + B*a*c + C*a*b)" and neuper@42376: equival_trans_4th_order: "(n/(a*b*c*d) = A/a + B/b + C/c + D/d) = neuper@42386: (n = A*b*c*d + B*a*c*d + C*a*b*d + D*a*b*c)" and neuper@42386: (*version 2: not yet used, see partial_fractions.sml*) neuper@42387: multiply_2nd_order: "(n/x = A/a + B/b) = (a*b*n/x = A*b + B*a)" and neuper@42387: multiply_3rd_order: "(n/x = A/a + B/b + C/c) = (a*b*c*n/x = A*b*c + B*a*c + C*a*b)" and neuper@42387: multiply_4th_order: neuper@42387: "(n/x = A/a + B/b + C/c + D/d) = (a*b*c*d*n/x = A*b*c*d + B*a*c*d + C*a*b*d + D*a*b*c)" neuper@42387: neuper@42387: text {* Probably the optimal formalization woudl be ... neuper@42387: neuper@42386: multiply_2nd_order: "x = a*b ==> (n/x = A/a + B/b) = (a*b*n/x = A*b + B*a)" and neuper@42386: multiply_3rd_order: "x = a*b*c ==> neuper@42386: (n/x = A/a + B/b + C/c) = (a*b*c*n/x = A*b*c + B*a*c + C*a*b)" and neuper@42386: multiply_4th_order: "x = a*b*c*d ==> neuper@42386: (n/x = A/a + B/b + C/c + D/d) = (a*b*c*d*n/x = A*b*c*d + B*a*c*d + C*a*b*d + D*a*b*c)" jan@42353: neuper@42387: ... because it would allow to start the ansatz as follows neuper@42387: (1) 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))) = 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))) neuper@42387: (2) 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))) = AA / (z - 1 / 2) + BB / (z - -1 / 4) neuper@42387: (3) (z - 1 / 2) * (z - -1 / 4) * 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))) = neuper@42387: (z - 1 / 2) * (z - -1 / 4) * AA / (z - 1 / 2) + BB / (z - -1 / 4) neuper@42387: (4) 3 = A * (z - -1 / 4) + B * (z - 1 / 2) neuper@42387: neuper@42387: ... (1==>2) ansatz neuper@42387: (2==>3) multiply_* neuper@42387: (3==>4) norm_Rational neuper@42387: TODOs for this version ar in partial_fractions.sml "--- progr.vers.2: " neuper@42387: *} neuper@42387: jan@42353: ML {* jan@42353: val ansatz_rls = prep_rls( jan@42353: Rls {id = "ansatz_rls", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@42451: erls = Erls, srls = Erls, calc = [], errpatts = [], jan@42353: rules = neuper@42376: [Thm ("ansatz_2nd_order",num_str @{thm ansatz_2nd_order}), neuper@42376: Thm ("ansatz_3rd_order",num_str @{thm ansatz_3rd_order}) jan@42353: ], jan@42353: scr = EmptyScr}:rls); jan@42353: jan@42354: val equival_trans = prep_rls( jan@42354: Rls {id = "equival_trans", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@42451: erls = Erls, srls = Erls, calc = [], errpatts = [], jan@42354: rules = neuper@42376: [Thm ("equival_trans_2nd_order",num_str @{thm equival_trans_2nd_order}), neuper@42376: Thm ("equival_trans_3rd_order",num_str @{thm equival_trans_3rd_order}) jan@42354: ], jan@42354: scr = EmptyScr}:rls); neuper@42386: neuper@42386: val multiply_ansatz = prep_rls( neuper@42386: Rls {id = "multiply_ansatz", preconds = [], rew_ord = ("dummy_ord",dummy_ord), neuper@42386: erls = Erls, neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@42386: rules = neuper@42386: [Thm ("multiply_2nd_order",num_str @{thm multiply_2nd_order}) neuper@42386: ], neuper@42386: scr = EmptyScr}:rls); jan@42354: *} jan@42354: jan@42354: text {*store the rule set for math engine*} neuper@52125: setup {* KEStore_Elems.add_rlss neuper@52125: [("ansatz_rls", (Context.theory_name @{theory}, ansatz_rls)), neuper@52125: ("multiply_ansatz", (Context.theory_name @{theory}, multiply_ansatz)), neuper@52125: ("equival_trans", (Context.theory_name @{theory}, equival_trans))] *} jan@42344: neuper@42376: subsection {* Specification *} jan@42344: neuper@42376: consts neuper@42376: decomposedFunction :: "real => una" neuper@42376: neuper@42289: ML {* neuper@42386: check_guhs_unique := false; (*WN120307 REMOVE after editing*) neuper@42289: *} s1210629013@55359: setup {* KEStore_Elems.add_pbts s1210629013@55339: [(prep_pbt @{theory} "pbl_simp_rat_partfrac" [] e_pblID s1210629013@55339: (["partial_fraction", "rational", "simplification"], s1210629013@55339: [("#Given" ,["functionTerm t_t", "solveFor v_v"]), s1210629013@55339: (* TODO: call this sub-problem with appropriate functionTerm: s1210629013@55339: leading coefficient of the denominator is 1: to be checked here! and.. s1210629013@55339: ("#Where" ,["((get_numerator t_t) has_degree_in v_v) < s1210629013@55339: ((get_denominator t_t) has_degree_in v_v)"]), TODO*) s1210629013@55339: ("#Find" ,["decomposedFunction p_p'''"])], s1210629013@55339: append_rls "e_rls" e_rls [(*for preds in where_ TODO*)], s1210629013@55339: NONE, s1210629013@55339: [["simplification","of_rationals","to_partial_fraction"]]))] *} jan@42354: neuper@42376: subsection {* Method *} neuper@42376: consts neuper@42376: PartFracScript :: "[real,real, real] => real" neuper@42376: ("((Script PartFracScript (_ _ =))// (_))" 9) jan@42353: neuper@42376: text {* rule set for functions called in the Script *} neuper@42376: ML {* neuper@42413: val srls_partial_fraction = Rls {id="srls_partial_fraction", neuper@42376: preconds = [], neuper@42376: rew_ord = ("termlessI",termlessI), neuper@42376: erls = append_rls "erls_in_srls_partial_fraction" e_rls neuper@42376: [(*for asm in NTH_CONS ...*) neuper@42376: Calc ("Orderings.ord_class.less",eval_equ "#less_"), neuper@42376: (*2nd NTH_CONS pushes n+-1 into asms*) neuper@42376: Calc("Groups.plus_class.plus", eval_binop "#add_")], neuper@42451: srls = Erls, calc = [], errpatts = [], neuper@42376: rules = [ neuper@42376: Thm ("NTH_CONS",num_str @{thm NTH_CONS}), neuper@42376: Calc("Groups.plus_class.plus", eval_binop "#add_"), neuper@42376: Thm ("NTH_NIL",num_str @{thm NTH_NIL}), neuper@42376: Calc("Tools.lhs", eval_lhs "eval_lhs_"), neuper@42376: Calc("Tools.rhs", eval_rhs"eval_rhs_"), neuper@42376: Calc("Atools.argument'_in", eval_argument_in "Atools.argument'_in"), neuper@42376: Calc("Rational.get_denominator", eval_get_denominator "#get_denominator"), neuper@42376: Calc("Rational.get_numerator", eval_get_numerator "#get_numerator"), neuper@42376: Calc("Partial_Fractions.factors_from_solution", neuper@42376: eval_factors_from_solution "#factors_from_solution"), neuper@42376: Calc("Partial_Fractions.drop_questionmarks", eval_drop_questionmarks "#drop_?")], neuper@42376: scr = EmptyScr}; neuper@42376: *} neuper@42376: ML {* neuper@42376: eval_drop_questionmarks; neuper@42376: *} neuper@42376: ML {* neuper@48761: val ctxt = Proof_Context.init_global @{theory}; neuper@42376: val SOME t = parseNEW ctxt "eqr = drop_questionmarks eqr"; neuper@42376: *} neuper@42376: ML {* neuper@42417: parseNEW ctxt "decomposedFunction p_p'''"; neuper@42376: parseNEW ctxt "decomposedFunction"; neuper@42376: *} neuper@42415: s1210629013@55380: (* current version, error outcommented *) s1210629013@55373: setup {* KEStore_Elems.add_mets s1210629013@55373: [prep_met @{theory} "met_partial_fraction" [] e_metID s1210629013@55373: (["simplification","of_rationals","to_partial_fraction"], s1210629013@55373: [("#Given" ,["functionTerm t_t", "solveFor v_v"]), s1210629013@55373: (*("#Where" ,["((get_numerator t_t) has_degree_in v_v) < s1210629013@55373: ((get_denominator t_t) has_degree_in v_v)"]), TODO*) s1210629013@55373: ("#Find" ,["decomposedFunction p_p'''"])], s1210629013@55373: (*f_f = 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z)), zzz: z*) s1210629013@55373: {rew_ord'="tless_true", rls'= e_rls, calc = [], srls = srls_partial_fraction, prls = e_rls, s1210629013@55373: crls = e_rls, errpats = [], nrls = e_rls}, s1210629013@55373: (*([], Frm), Problem (Partial_Fractions, [partial_fraction, rational, simplification])*) s1210629013@55373: "Script PartFracScript (f_f::real) (zzz::real) = " ^ s1210629013@55373: (*([1], Frm), 3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z)))*) s1210629013@55373: "(let f_f = Take f_f; " ^ s1210629013@55373: (* num_orig = 3*) s1210629013@55373: " (num_orig::real) = get_numerator f_f; " ^ s1210629013@55373: (*([1], Res), 24 / (-1 + -2 * z + 8 * z ^^^ 2)*) s1210629013@55373: " f_f = (Rewrite_Set norm_Rational False) f_f; " ^ s1210629013@55373: (* denom = -1 + -2 * z + 8 * z ^^^ 2*) s1210629013@55373: " (denom::real) = get_denominator f_f; " ^ s1210629013@55373: (* equ = -1 + -2 * z + 8 * z ^^^ 2 = 0*) s1210629013@55373: " (equ::bool) = (denom = (0::real)); " ^ s1210629013@55373: s1210629013@55373: (*([2], Pbl), solve (-1 + -2 * z + 8 * z ^^^ 2 = 0, z)*) s1210629013@55373: " (L_L::bool list) = (SubProblem (PolyEq', " ^ s1210629013@55373: " [abcFormula, degree_2, polynomial, univariate, equation], " ^ s1210629013@55373: (*([2], Res), [z = 1 / 2, z = -1 / 4]*) s1210629013@55373: " [no_met]) [BOOL equ, REAL zzz]); " ^ s1210629013@55373: (* facs: (z - 1 / 2) * (z - -1 / 4)*) s1210629013@55373: " (facs::real) = factors_from_solution L_L; " ^ s1210629013@55373: (*([3], Frm), 33 / ((z - 1 / 2) * (z - -1 / 4)) *) s1210629013@55373: " (eql::real) = Take (num_orig / facs); " ^ s1210629013@55373: (*([3], Res), ?A / (z - 1 / 2) + ?B / (z - -1 / 4)*) s1210629013@55373: " (eqr::real) = (Try (Rewrite_Set ansatz_rls False)) eql; " ^ s1210629013@55373: (*([4], Frm), 3 / ((z - 1 / 2) * (z - -1 / 4)) = ?A / (z - 1 / 2) + ?B / (z - -1 / 4)*) s1210629013@55373: " (eq::bool) = Take (eql = eqr); " ^ s1210629013@55373: (*([4], Res), 3 = ?A * (z - -1 / 4) + ?B * (z - 1 / 2)*) s1210629013@55373: " eq = (Try (Rewrite_Set equival_trans False)) eq;" ^ s1210629013@55373: (* eq = 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*) s1210629013@55373: " eq = drop_questionmarks eq; " ^ s1210629013@55373: (* z1 = 1 / 2*) s1210629013@55373: " (z1::real) = (rhs (NTH 1 L_L)); " ^ s1210629013@55373: (* z2 = -1 / 4*) s1210629013@55373: " (z2::real) = (rhs (NTH 2 L_L)); " ^ s1210629013@55373: (*([5], Frm), 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*) s1210629013@55373: " (eq_a::bool) = Take eq; " ^ s1210629013@55373: (*([5], Res), 3 = A * (1 / 2 - -1 / 4) + B * (1 / 2 - 1 / 2)*) s1210629013@55373: " eq_a = (Substitute [zzz = z1]) eq; " ^ s1210629013@55373: (*([6], Res), 3 = 3 * A / 4*) s1210629013@55373: " eq_a = (Rewrite_Set norm_Rational False) eq_a; " ^ s1210629013@55373: s1210629013@55373: (*([7], Pbl), solve (3 = 3 * A / 4, A)*) s1210629013@55373: " (sol_a::bool list) = " ^ s1210629013@55373: " (SubProblem (Isac', [univariate,equation], [no_met]) " ^ s1210629013@55373: (*([7], Res), [A = 4]*) s1210629013@55373: " [BOOL eq_a, REAL (A::real)]); " ^ s1210629013@55373: (* a = 4*) s1210629013@55373: " (a::real) = (rhs (NTH 1 sol_a)); " ^ s1210629013@55373: (*([8], Frm), 3 = A * (z - -1 / 4) + B * (z - 1 / 2)*) s1210629013@55373: " (eq_b::bool) = Take eq; " ^ s1210629013@55373: (*([8], Res), 3 = A * (-1 / 4 - -1 / 4) + B * (-1 / 4 - 1 / 2)*) s1210629013@55373: " eq_b = (Substitute [zzz = z2]) eq_b; " ^ s1210629013@55373: (*([9], Res), 3 = -3 * B / 4*) s1210629013@55373: " eq_b = (Rewrite_Set norm_Rational False) eq_b; " ^ s1210629013@55373: (*([10], Pbl), solve (3 = -3 * B / 4, B)*) s1210629013@55373: " (sol_b::bool list) = " ^ s1210629013@55373: " (SubProblem (Isac', [univariate,equation], [no_met]) " ^ s1210629013@55373: (*([10], Res), [B = -4]*) s1210629013@55373: " [BOOL eq_b, REAL (B::real)]); " ^ s1210629013@55373: (* b = -4*) s1210629013@55373: " (b::real) = (rhs (NTH 1 sol_b)); " ^ s1210629013@55373: (* eqr = A / (z - 1 / 2) + B / (z - -1 / 4)*) s1210629013@55373: " eqr = drop_questionmarks eqr; " ^ s1210629013@55373: (*([11], Frm), A / (z - 1 / 2) + B / (z - -1 / 4)*) s1210629013@55373: " (pbz::real) = Take eqr; " ^ s1210629013@55373: (*([11], Res), 4 / (z - 1 / 2) + -4 / (z - -1 / 4)*) s1210629013@55373: " pbz = ((Substitute [A = a, B = b]) pbz) " ^ s1210629013@55373: (*([], Res), 4 / (z - 1 / 2) + -4 / (z - -1 / 4)*) s1210629013@55373: "in pbz)" s1210629013@55373: )] s1210629013@55373: *} neuper@42376: ML {* neuper@42376: (* neuper@42376: val fmz = neuper@42376: ["functionTerm (3 / (z * (z - 1 / 4 + -1 / 8 * (1 / z))))", neuper@42376: "solveFor z", "functionTerm p_p"]; neuper@42376: val (dI',pI',mI') = neuper@42376: ("Partial_Fractions", neuper@42376: ["partial_fraction", "rational", "simplification"], neuper@42376: ["simplification","of_rationals","to_partial_fraction"]); neuper@42376: val (p,_,f,nxt,_,pt) = CalcTreeTEST [(fmz, (dI',pI',mI'))]; neuper@42376: *) neuper@42376: *} neuper@42289: jan@42295: neuper@42376: neuper@42376: subsection {**} neuper@42376: neuper@42289: end