neuper@37906: (* tools for integration over the reals
neuper@37906:    author: Walther Neuper 050905, 08:51
neuper@37906:    (c) due to copyright terms
neuper@37906: 
neuper@37947: use"Knowledge/Integrate.ML";
neuper@37906: use"Integrate.ML";
neuper@37906: 
neuper@37906: remove_thy"Integrate";
neuper@37947: use_thy"Knowledge/Isac";
neuper@37906: *)
neuper@37906: 
neuper@37906: (** interface isabelle -- isac **)
neuper@37906: 
neuper@37906: theory' := overwritel (!theory', [("Integrate.thy",Integrate.thy)]);
neuper@37906: 
neuper@37906: (** eval functions **)
neuper@37906: 
neuper@37906: val c = Free ("c", HOLogic.realT);
neuper@37906: (*.create a new unique variable 'c..' in a term; for use by Calc in a rls;
neuper@37906:    an alternative to do this would be '(Try (Calculate new_c_) (new_c es__))'
neuper@37906:    in the script; this will be possible if currying doesnt take the value
neuper@37906:    from a variable, but the value '(new_c es__)' itself.*)
neuper@37906: fun new_c term = 
neuper@37906:     let fun selc var = 
neuper@37906: 	    case (explode o id_of) var of
neuper@37906: 		"c"::[] => true
neuper@37906: 	      |	"c"::"_"::is => (case (int_of_str o implode) is of
neuper@37926: 				     SOME _ => true
neuper@37926: 				   | NONE => false)
neuper@37906:               | _ => false;
neuper@37906: 	fun get_coeff c = case (explode o id_of) c of
neuper@37906: 	      		      "c"::"_"::is => (the o int_of_str o implode) is
neuper@37906: 			    | _ => 0;
neuper@37906:         val cs = filter selc (vars term);
neuper@37906:     in 
neuper@37906: 	case cs of
neuper@37906: 	    [] => c
neuper@37906: 	  | [c] => Free ("c_2", HOLogic.realT)
neuper@37906: 	  | cs => 
neuper@37906: 	    let val max_coeff = maxl (map get_coeff cs)
neuper@37906: 	    in Free ("c_"^string_of_int (max_coeff + 1), HOLogic.realT) end
neuper@37906:     end;
neuper@37906: 
neuper@37906: (*WN080222
neuper@37906: (*("new_c", ("Integrate.new'_c", eval_new_c "#new_c_"))*)
neuper@37906: fun eval_new_c _ _ (p as (Const ("Integrate.new'_c",_) $ t)) _ =
neuper@37926:      SOME ((term2str p) ^ " = " ^ term2str (new_c p),
neuper@37906: 	  Trueprop $ (mk_equality (p, new_c p)))
neuper@37926:   | eval_new_c _ _ _ _ = NONE;
neuper@37906: *)
neuper@37906: 
neuper@37906: (*WN080222:*)
neuper@37906: (*("add_new_c", ("Integrate.add'_new'_c", eval_add_new_c "#add_new_c_"))
neuper@37906:   add a new c to a term or a fun-equation;
neuper@37906:   this is _not in_ the term, because only applied to _whole_ term*)
neuper@37906: fun eval_add_new_c (_:string) "Integrate.add'_new'_c" p (_:theory) =
neuper@37906:     let val p' = case p of
neuper@37906: 		     Const ("op =", T) $ lh $ rh => 
neuper@37906: 		     Const ("op =", T) $ lh $ mk_add rh (new_c rh)
neuper@37906: 		   | p => mk_add p (new_c p)
neuper@37926:     in SOME ((term2str p) ^ " = " ^ term2str p',
neuper@37906: 	  Trueprop $ (mk_equality (p, p')))
neuper@37906:     end
neuper@37926:   | eval_add_new_c _ _ _ _ = NONE;
neuper@37906: 
neuper@37906: 
neuper@37906: (*("is_f_x", ("Integrate.is'_f'_x", eval_is_f_x "is_f_x_"))*)
neuper@37906: fun eval_is_f_x _ _(p as (Const ("Integrate.is'_f'_x", _)
neuper@37906: 					   $ arg)) _ =
neuper@37906:     if is_f_x arg
neuper@37926:     then SOME ((term2str p) ^ " = True",
neuper@37906: 	       Trueprop $ (mk_equality (p, HOLogic.true_const)))
neuper@37926:     else SOME ((term2str p) ^ " = False",
neuper@37906: 	       Trueprop $ (mk_equality (p, HOLogic.false_const)))
neuper@37926:   | eval_is_f_x _ _ _ _ = NONE;
neuper@37906: 
neuper@37906: calclist':= overwritel (!calclist', 
neuper@37906:    [(*("new_c", ("Integrate.new'_c", eval_new_c "new_c_")),*)
neuper@37906:     ("add_new_c", ("Integrate.add'_new'_c", eval_add_new_c "add_new_c_")),
neuper@37906:     ("is_f_x", ("Integrate.is'_f'_x", eval_is_f_x "is_f_idextifier_"))
neuper@37906:     ]);
neuper@37906: 
neuper@37906: 
neuper@37906: (** rulesets **)
neuper@37906: 
neuper@37906: (*.rulesets for integration.*)
neuper@37906: val integration_rules = 
neuper@37906:     Rls {id="integration_rules", preconds = [], 
neuper@37906: 	 rew_ord = ("termlessI",termlessI), 
neuper@37906: 	 erls = Rls {id="conditions_in_integration_rules", 
neuper@37906: 		     preconds = [], 
neuper@37906: 		     rew_ord = ("termlessI",termlessI), 
neuper@37906: 		     erls = Erls, 
neuper@37906: 		     srls = Erls, calc = [],
neuper@37906: 		     rules = [(*for rewriting conditions in Thm's*)
neuper@37906: 			      Calc ("Atools.occurs'_in", 
neuper@37906: 				    eval_occurs_in "#occurs_in_"),
neuper@37906: 			      Thm ("not_true",num_str not_true),
neuper@37906: 			      Thm ("not_false",not_false)
neuper@37906: 			      ],
neuper@37906: 		     scr = EmptyScr}, 
neuper@37906: 	 srls = Erls, calc = [],
neuper@37906: 	 rules = [
neuper@37906: 		  Thm ("integral_const",num_str integral_const),
neuper@37906: 		  Thm ("integral_var",num_str integral_var),
neuper@37906: 		  Thm ("integral_add",num_str integral_add),
neuper@37906: 		  Thm ("integral_mult",num_str integral_mult),
neuper@37906: 		  Thm ("integral_pow",num_str integral_pow),
neuper@37906: 		  Calc ("op +", eval_binop "#add_")(*for n+1*)
neuper@37906: 		  ],
neuper@37906: 	 scr = EmptyScr};
neuper@37906: val add_new_c = 
neuper@37906:     Seq {id="add_new_c", preconds = [], 
neuper@37906: 	 rew_ord = ("termlessI",termlessI), 
neuper@37906: 	 erls = Rls {id="conditions_in_add_new_c", 
neuper@37906: 		     preconds = [], 
neuper@37906: 		     rew_ord = ("termlessI",termlessI), 
neuper@37906: 		     erls = Erls, 
neuper@37906: 		     srls = Erls, calc = [],
neuper@37906: 		     rules = [Calc ("Tools.matches", eval_matches""),
neuper@37906: 			      Calc ("Integrate.is'_f'_x", 
neuper@37906: 				    eval_is_f_x "is_f_x_"),
neuper@37906: 			      Thm ("not_true",num_str not_true),
neuper@37906: 			      Thm ("not_false",num_str not_false)
neuper@37906: 			      ],
neuper@37906: 		     scr = EmptyScr}, 
neuper@37906: 	 srls = Erls, calc = [],
neuper@37906: 	 rules = [ (*Thm ("call_for_new_c", num_str call_for_new_c),*)
neuper@37906: 		   Cal1 ("Integrate.add'_new'_c", eval_add_new_c "new_c_")
neuper@37906: 		   ],
neuper@37906: 	 scr = EmptyScr};
neuper@37906: 
neuper@37906: (*.rulesets for simplifying Integrals.*)
neuper@37906: 
neuper@37906: (*.for simplify_Integral adapted from 'norm_Rational_rls'.*)
neuper@37906: val norm_Rational_rls_noadd_fractions = 
neuper@37906: Rls {id = "norm_Rational_rls_noadd_fractions", preconds = [], 
neuper@37906:      rew_ord = ("dummy_ord",dummy_ord), 
neuper@37906:      erls = norm_rat_erls, srls = Erls, calc = [],
neuper@37906:      rules = [(*Rls_ common_nominator_p_rls,!!!*)
neuper@37906: 	      Rls_ (*rat_mult_div_pow original corrected WN051028*)
neuper@37906: 		  (Rls {id = "rat_mult_div_pow", preconds = [], 
neuper@37906: 		       rew_ord = ("dummy_ord",dummy_ord), 
neuper@37906: 		       erls = (*FIXME.WN051028 e_rls,*)
neuper@37906: 		       append_rls "e_rls-is_polyexp" e_rls
neuper@37906: 				  [Calc ("Poly.is'_polyexp", 
neuper@37906: 					 eval_is_polyexp "")],
neuper@37906: 				  srls = Erls, calc = [],
neuper@37906: 				  rules = [Thm ("rat_mult",num_str rat_mult),
neuper@37906: 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)
neuper@37906: 	       Thm ("rat_mult_poly_l",num_str rat_mult_poly_l),
neuper@37906: 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)
neuper@37906: 	       Thm ("rat_mult_poly_r",num_str rat_mult_poly_r),
neuper@37906: 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)
neuper@37906: 
neuper@37906: 	       Thm ("real_divide_divide1_mg", real_divide_divide1_mg),
neuper@37906: 	       (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)
neuper@37906: 	       Thm ("real_divide_divide1_eq", real_divide_divide1_eq),
neuper@37906: 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)
neuper@37906: 	       Thm ("real_divide_divide2_eq", real_divide_divide2_eq),
neuper@37906: 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)
neuper@37906: 	       Calc ("HOL.divide"  ,eval_cancel "#divide_"),
neuper@37906: 	      
neuper@37906: 	       Thm ("rat_power", num_str rat_power)
neuper@37906: 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)
neuper@37906: 	       ],
neuper@37906:       scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37906:       }),
neuper@37906: 		Rls_ make_rat_poly_with_parentheses,
neuper@37906: 		Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)
neuper@37906: 		Rls_ rat_reduce_1
neuper@37906: 		],
neuper@37906:        scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37906:        }:rls;
neuper@37906: 
neuper@37906: (*.for simplify_Integral adapted from 'norm_Rational'.*)
neuper@37906: val norm_Rational_noadd_fractions = 
neuper@37906:    Seq {id = "norm_Rational_noadd_fractions", preconds = [], 
neuper@37906:        rew_ord = ("dummy_ord",dummy_ord), 
neuper@37906:        erls = norm_rat_erls, srls = Erls, calc = [],
neuper@37906:        rules = [Rls_ discard_minus_,
neuper@37906: 		Rls_ rat_mult_poly,(* removes double fractions like a/b/c    *)
neuper@37906: 		Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)
neuper@37906: 		Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)
neuper@37906: 		Rls_ norm_Rational_rls_noadd_fractions,(* the main rls (#)   *)
neuper@37906: 		Rls_ discard_parentheses_ (* mult only                       *)
neuper@37906: 		],
neuper@37906:        scr = Script ((term_of o the o (parse thy)) "empty_script")
neuper@37906:        }:rls;
neuper@37906: 
neuper@37906: (*.simplify terms before and after Integration such that  
neuper@37906:    ..a.x^2/2 + b.x^3/3.. is made to ..a/2.x^2 + b/3.x^3.. (and NO
neuper@37906:    common denominator as done by norm_Rational or make_ratpoly_in.
neuper@37906:    This is a copy from 'make_ratpoly_in' with respective reduction of rules and
neuper@37906:    *1* expand the term, ie. distribute * and / over +
neuper@37906: .*)
neuper@37906: val separate_bdv2 =
neuper@37906:     append_rls "separate_bdv2"
neuper@37906: 	       collect_bdv
neuper@37906: 	       [Thm ("separate_bdv", num_str separate_bdv),
neuper@37906: 		(*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*)
neuper@37906: 		Thm ("separate_bdv_n", num_str separate_bdv_n),
neuper@37906: 		Thm ("separate_1_bdv", num_str separate_1_bdv),
neuper@37906: 		(*"?bdv / ?b = (1 / ?b) * ?bdv"*)
neuper@37906: 		Thm ("separate_1_bdv_n", num_str separate_1_bdv_n)(*,
neuper@37906: 			  (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*)
neuper@37906: 			  *****Thm ("real_add_divide_distrib", 
neuper@37906: 			  *****num_str real_add_divide_distrib)
neuper@37906: 			  (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)----------*)
neuper@37906: 		];
neuper@37906: val simplify_Integral = 
neuper@37906:   Seq {id = "simplify_Integral", preconds = []:term list, 
neuper@37906:        rew_ord = ("dummy_ord", dummy_ord),
neuper@37906:       erls = Atools_erls, srls = Erls,
neuper@37906:       calc = [], (*asm_thm = [],*)
neuper@37906:       rules = [Thm ("real_add_mult_distrib",num_str real_add_mult_distrib),
neuper@37906:  	       (*"(?z1.0 + ?z2.0) * ?w = ?z1.0 * ?w + ?z2.0 * ?w"*)
neuper@37906: 	       Thm ("real_add_divide_distrib",num_str real_add_divide_distrib),
neuper@37906:  	       (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)
neuper@37906: 	       (*^^^^^ *1* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*)
neuper@37906: 	       Rls_ norm_Rational_noadd_fractions,
neuper@37906: 	       Rls_ order_add_mult_in,
neuper@37906: 	       Rls_ discard_parentheses,
neuper@37906: 	       (*Rls_ collect_bdv, from make_polynomial_in*)
neuper@37906: 	       Rls_ separate_bdv2,
neuper@37906: 	       Calc ("HOL.divide"  ,eval_cancel "#divide_")
neuper@37906: 	       ],
neuper@37906:       scr = EmptyScr}:rls;      
neuper@37906: 
neuper@37906: 
neuper@37906: (*simplify terms before and after Integration such that  
neuper@37906:    ..a.x^2/2 + b.x^3/3.. is made to ..a/2.x^2 + b/3.x^3.. (and NO
neuper@37906:    common denominator as done by norm_Rational or make_ratpoly_in.
neuper@37906:    This is a copy from 'make_polynomial_in' with insertions from 
neuper@37906:    'make_ratpoly_in' 
neuper@37906: THIS IS KEPT FOR COMPARISON ............................................   
neuper@37906: * val simplify_Integral = prep_rls(
neuper@37906: *   Seq {id = "", preconds = []:term list, 
neuper@37906: *        rew_ord = ("dummy_ord", dummy_ord),
neuper@37906: *       erls = Atools_erls, srls = Erls,
neuper@37906: *       calc = [], (*asm_thm = [],*)
neuper@37906: *       rules = [Rls_ expand_poly,
neuper@37906: * 	       Rls_ order_add_mult_in,
neuper@37906: * 	       Rls_ simplify_power,
neuper@37906: * 	       Rls_ collect_numerals,
neuper@37906: * 	       Rls_ reduce_012,
neuper@37906: * 	       Thm ("realpow_oneI",num_str realpow_oneI),
neuper@37906: * 	       Rls_ discard_parentheses,
neuper@37906: * 	       Rls_ collect_bdv,
neuper@37906: * 	       (*below inserted from 'make_ratpoly_in'*)
neuper@37906: * 	       Rls_ (append_rls "separate_bdv"
neuper@37906: * 			 collect_bdv
neuper@37906: * 			 [Thm ("separate_bdv", num_str separate_bdv),
neuper@37906: * 			  (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*)
neuper@37906: * 			  Thm ("separate_bdv_n", num_str separate_bdv_n),
neuper@37906: * 			  Thm ("separate_1_bdv", num_str separate_1_bdv),
neuper@37906: * 			  (*"?bdv / ?b = (1 / ?b) * ?bdv"*)
neuper@37906: * 			  Thm ("separate_1_bdv_n", num_str separate_1_bdv_n)(*,
neuper@37906: * 			  (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*)
neuper@37906: * 			  Thm ("real_add_divide_distrib", 
neuper@37906: * 				 num_str real_add_divide_distrib)
neuper@37906: * 			   (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)*)
neuper@37906: * 			  ]),
neuper@37906: * 	       Calc ("HOL.divide"  ,eval_cancel "#divide_")
neuper@37906: * 	       ],
neuper@37906: *       scr = EmptyScr
neuper@37906: *       }:rls); 
neuper@37906: .......................................................................*)
neuper@37906: 
neuper@37906: val integration = 
neuper@37906:     Seq {id="integration", preconds = [], 
neuper@37906: 	 rew_ord = ("termlessI",termlessI), 
neuper@37906: 	 erls = Rls {id="conditions_in_integration", 
neuper@37906: 		     preconds = [], 
neuper@37906: 		     rew_ord = ("termlessI",termlessI), 
neuper@37906: 		     erls = Erls, 
neuper@37906: 		     srls = Erls, calc = [],
neuper@37906: 		     rules = [],
neuper@37906: 		     scr = EmptyScr}, 
neuper@37906: 	 srls = Erls, calc = [],
neuper@37906: 	 rules = [ Rls_ integration_rules,
neuper@37906: 		   Rls_ add_new_c,
neuper@37906: 		   Rls_ simplify_Integral
neuper@37906: 		   ],
neuper@37906: 	 scr = EmptyScr};
neuper@37906: ruleset' := 
neuper@37906: overwritelthy thy (!ruleset', 
neuper@37906: 	    [("integration_rules", prep_rls integration_rules),
neuper@37906: 	     ("add_new_c", prep_rls add_new_c),
neuper@37906: 	     ("simplify_Integral", prep_rls simplify_Integral),
neuper@37906: 	     ("integration", prep_rls integration),
neuper@37906: 	     ("separate_bdv2", separate_bdv2),
neuper@37906: 	     ("norm_Rational_noadd_fractions", norm_Rational_noadd_fractions),
neuper@37906: 	     ("norm_Rational_rls_noadd_fractions", 
neuper@37906: 	      norm_Rational_rls_noadd_fractions)
neuper@37906: 	     ]);
neuper@37906: 
neuper@37906: (** problems **)
neuper@37906: 
neuper@37906: store_pbt
neuper@37906:  (prep_pbt Integrate.thy "pbl_fun_integ" [] e_pblID
neuper@37906:  (["integrate","function"],
neuper@37906:   [("#Given" ,["functionTerm f_", "integrateBy v_"]),
neuper@37906:    ("#Find"  ,["antiDerivative F_"])
neuper@37906:   ],
neuper@37906:   append_rls "e_rls" e_rls [(*for preds in where_*)], 
neuper@37926:   SOME "Integrate (f_, v_)", 
neuper@37906:   [["diff","integration"]]));
neuper@37906:  
neuper@37906: (*here "named" is used differently from Differentiation"*)
neuper@37906: store_pbt
neuper@37906:  (prep_pbt Integrate.thy "pbl_fun_integ_nam" [] e_pblID
neuper@37906:  (["named","integrate","function"],
neuper@37906:   [("#Given" ,["functionTerm f_", "integrateBy v_"]),
neuper@37906:    ("#Find"  ,["antiDerivativeName F_"])
neuper@37906:   ],
neuper@37906:   append_rls "e_rls" e_rls [(*for preds in where_*)], 
neuper@37926:   SOME "Integrate (f_, v_)", 
neuper@37906:   [["diff","integration","named"]]));
neuper@37906:  
neuper@37906: (** methods **)
neuper@37906: 
neuper@37906: store_met
neuper@37906:     (prep_met Integrate.thy "met_diffint" [] e_metID
neuper@37906: 	      (["diff","integration"],
neuper@37906: 	       [("#Given" ,["functionTerm f_", "integrateBy v_"]),
neuper@37906: 		("#Find"  ,["antiDerivative F_"])
neuper@37906: 		],
neuper@37906: 	       {rew_ord'="tless_true", rls'=Atools_erls, calc = [], 
neuper@37906: 		srls = e_rls, 
neuper@37906: 		prls=e_rls,
neuper@37906: 	     crls = Atools_erls, nrls = e_rls},
neuper@37906: "Script IntegrationScript (f_::real) (v_::real) =                \
neuper@37906: \  (let t_ = Take (Integral f_ D v_)                             \
neuper@37906: \   in (Rewrite_Set_Inst [(bdv,v_)] integration False) (t_::real))"
neuper@37906: ));
neuper@37906:     
neuper@37906: store_met
neuper@37906:     (prep_met Integrate.thy "met_diffint_named" [] e_metID
neuper@37906: 	      (["diff","integration","named"],
neuper@37906: 	       [("#Given" ,["functionTerm f_", "integrateBy v_"]),
neuper@37906: 		("#Find"  ,["antiDerivativeName F_"])
neuper@37906: 		],
neuper@37906: 	       {rew_ord'="tless_true", rls'=Atools_erls, calc = [], 
neuper@37906: 		srls = e_rls, 
neuper@37906: 		prls=e_rls,
neuper@37906: 		crls = Atools_erls, nrls = e_rls},
neuper@37906: "Script NamedIntegrationScript (f_::real) (v_::real) (F_::real=>real) = \
neuper@37906: \  (let t_ = Take (F_ v_ = Integral f_ D v_)                         \
neuper@37906: \   in ((Try (Rewrite_Set_Inst [(bdv,v_)] simplify_Integral False)) @@\
neuper@37906: \       (Rewrite_Set_Inst [(bdv,v_)] integration False)) t_)"
neuper@37906: (*
neuper@37906: "Script NamedIntegrationScript (f_::real) (v_::real) (F_::real=>real) = \
neuper@37906: \  (let t_ = Take (F_ v_ = Integral f_ D v_)                         \
neuper@37906: \   in (Rewrite_Set_Inst [(bdv,v_)] integration False) t_)"
neuper@37906: *)
neuper@37906:  ));