(* integration over the reals

    author: Walther Neuper

    050814, 08:51

    (c) due to copyright terms

 *)

 theory Integrate imports Diff begin

 consts

   Integral            :: "[real, real]=> real" ("Integral _ D _" 91)

   add_new_c          :: "real => real"        ("add'_new'_c _" 66) 

   is_f_x            :: "real => bool"        ("_ is'_f'_x" 10)

   (*descriptions in the related problems*)

   integrateBy         :: "real => una"

   antiDerivative      :: "real => una"

   antiDerivativeName  :: "(real => real) => una"

   (*the CAS-command, eg. "Integrate (2*x \<up> 3, x)"*)

   Integrate           :: "[real * real] => real"

 axiomatization where

 (*stated as axioms, todo: prove as theorems

   'bdv' is a constant handled on the meta-level 

    specifically as a 'bound variable'            *)

 (*Ambiguous input\<^here> produces 3 parse trees -----------------------------\\*)

   integral_const:    "Not (bdv occurs_in u) ==> Integral u D bdv = u * bdv" and

   integral_var:      "Integral bdv D bdv = bdv \<up> 2 / 2" and

   integral_add:      "Integral (u + v) D bdv =  

 		     (Integral u D bdv) + (Integral v D bdv)" and

   integral_mult:     "[| Not (bdv occurs_in u); bdv occurs_in v |] ==>  

 		     Integral (u * v) D bdv = u * (Integral v D bdv)" and

 (*WN080222: this goes into sub-terms, too ...

   call_for_new_c:    "[| Not (matches (u + new_c v) a); Not (a is_f_x) |] ==>  

 		     a = a + new_c a"

 *)

   integral_pow:      "Integral bdv \<up> n D bdv = bdv \<up> (n+1) / (n + 1)"

 (*Ambiguous input\<^here> produces 3 parse trees -----------------------------//*)

 ML \<open>

 (** eval functions **)

 val c = Free ("c", HOLogic.realT);

 (*.create a new unique variable 'c..' in a term; for use by \<^rule_eval> in a rls;

    an alternative to do this would be '(Try (Calculate new_c_) (new_c es__))'

    in the program; this would be possible if currying would not take the value

    from a variable, but the value '(new_c es__)' itself.*)

 fun new_c term = 

   let

     fun selc var = 

 	    case (Symbol.explode o TermC.id_of) var of

 		    "c"::[] => true

 	    |	"c" :: "_":: is =>

         (case (TermC.int_opt_of_string o implode) is of

 				   SOME _ => true

 				 | NONE => false)

          | _ => false;

 	  fun get_coeff c =

       case (Symbol.explode o TermC.id_of) c of

 	      "c"::"_"::is => (the o TermC.int_opt_of_string o implode) is

 			| _ => 0;

     val cs = filter selc (TermC.vars term);

   in 

 	  case cs of

 	    [] => c

 	  | [_] => Free ("c_2", HOLogic.realT)

 	  | cs => 

 	    let val max_coeff = maxl (map get_coeff cs)

 	    in Free ("c_"^string_of_int (max_coeff + 1), HOLogic.realT) end

   end;

 (*WN080222:*)

 (*("add_new_c", ("Integrate.add_new_c", eval_add_new_c "#add_new_c_"))

   add a new c to a term or a fun-equation;

   this is _not in_ the term, because only applied to _whole_ term*)

 fun eval_add_new_c (_:string) "Integrate.add_new_c" p (_:Proof.context) =

     let val p' = case p of

 		     Const (\<^const_name>\<open>HOL.eq\<close>, T) $ lh $ rh => 

 		     Const (\<^const_name>\<open>HOL.eq\<close>, T) $ lh $ TermC.mk_add rh (new_c rh)

 		   | p => TermC.mk_add p (new_c p)

     in SOME ((UnparseC.term @{context} p) ^ " = " ^ UnparseC.term @{context} p',

 	  HOLogic.Trueprop $ (TermC.mk_equality (p, p')))

     end

   | eval_add_new_c _ _ _ _ = NONE;

 (*("is_f_x", ("Integrate.is_f_x", eval_is_f_x "is_f_x_"))*)

 fun eval_is_f_x _ _(p as (Const (\<^const_name>\<open>Integrate.is_f_x\<close>, _)

 					   $ arg)) _ =

     if TermC.is_f_x arg

     then SOME ((UnparseC.term @{context} p) ^ " = True",

 	       HOLogic.Trueprop $ (TermC.mk_equality (p, @{term True})))

     else SOME ((UnparseC.term @{context} p) ^ " = False",

 	       HOLogic.Trueprop $ (TermC.mk_equality (p, @{term False})))

   | eval_is_f_x _ _ _ _ = NONE;

 \<close>

 calculation add_new_c = \<open>eval_add_new_c "add_new_c_"\<close>

 calculation is_f_x = \<open>eval_is_f_x "is_f_idextifier_"\<close>

 ML \<open>

 (** rulesets **)

 (*.rulesets for integration.*)

 val integration_rules = 

   Rule_Def.Repeat {id="integration_rules", preconds = [], 

     rew_ord = ("termlessI",termlessI), 

     asm_rls = Rule_Def.Repeat {id="conditions_in_integration_rules", 

    	  preconds = [], 

    	  rew_ord = ("termlessI",termlessI), 

    	  asm_rls = Rule_Set.Empty, 

    	  prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

    	  rules = [(*for rewriting conditions in Thm's*)

    		   \<^rule_eval>\<open>Prog_Expr.occurs_in\<close> (Prog_Expr.eval_occurs_in "#occurs_in_"),

    		   \<^rule_thm>\<open>not_true\<close>,

    		   \<^rule_thm>\<open>not_false\<close>],

    	  program = Rule.Empty_Prog}, 

     prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

     rules = [

    	  \<^rule_thm>\<open>integral_const\<close>,

    	  \<^rule_thm>\<open>integral_var\<close>,

    	  \<^rule_thm>\<open>integral_add\<close>,

    	  \<^rule_thm>\<open>integral_mult\<close>,

    	  \<^rule_thm>\<open>integral_pow\<close>,

    	  \<^rule_eval>\<open>plus\<close> (**)(Calc_Binop.numeric "#add_")(*for n+1*)],

     program = Rule.Empty_Prog};

 \<close>

 ML \<open>

 val add_new_c = 

   Rule_Set.Sequence {id="add_new_c", preconds = [], 

     rew_ord = ("termlessI",termlessI), 

     asm_rls = Rule_Def.Repeat {id="conditions_in_add_new_c", 

       preconds = [], rew_ord = ("termlessI",termlessI), asm_rls = Rule_Set.Empty, 

       prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

       rules = [

         \<^rule_eval>\<open>Prog_Expr.matches\<close> (Prog_Expr.eval_matches""),

    	    \<^rule_eval>\<open>Integrate.is_f_x\<close> (eval_is_f_x "is_f_x_"),

    	    \<^rule_thm>\<open>not_true\<close>,

    	    \<^rule_thm>\<open>not_false\<close>],

       program = Rule.Empty_Prog}, 

     prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

     rules = [ (*\<^rule_thm>\<open>call_for_new_c\<close>,*)

       Rule.Cal1 ("Integrate.add_new_c", eval_add_new_c "new_c_")],

     program = Rule.Empty_Prog};

 \<close>

 ML \<open>

 (*.rulesets for simplifying Integrals.*)

 (*.for simplify_Integral adapted from 'norm_Rational_rls'.*)

 val norm_Rational_rls_noadd_fractions = 

   Rule_Def.Repeat {id = "norm_Rational_rls_noadd_fractions", preconds = [], 

     rew_ord = ("dummy_ord",Rewrite_Ord.function_empty), 

     asm_rls = norm_rat_erls, prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

     rules = [(*Rule.Rls_ add_fractions_p_rls,!!!*)

   	  Rule.Rls_ (*rat_mult_div_pow original corrected WN051028*)

   		(Rule_Def.Repeat {id = "rat_mult_div_pow", preconds = [], 

   		   rew_ord = ("dummy_ord",Rewrite_Ord.function_empty), 

   		   asm_rls = Rule_Set.append_rules "Rule_Set.empty-is_polyexp" Rule_Set.empty

   				 [\<^rule_eval>\<open>is_polyexp\<close> (eval_is_polyexp "")],

   			 prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

   		   rules = [

            \<^rule_thm>\<open>rat_mult\<close>, (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

   	       \<^rule_thm>\<open>rat_mult_poly_l\<close>, (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

   	       \<^rule_thm>\<open>rat_mult_poly_r\<close>, (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

   	       \<^rule_thm>\<open>real_divide_divide1_mg\<close>, (*"y ~= 0 ==> (u / v) / (y / z) = (u * z) / (y * v)"*)

   	       \<^rule_thm>\<open>divide_divide_eq_right\<close>, (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

   	       \<^rule_thm>\<open>divide_divide_eq_left\<close>, (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

   	       \<^rule_eval>\<open>divide\<close> (Prog_Expr.eval_cancel "#divide_e"),

   	       \<^rule_thm>\<open>rat_power\<close>], (*"(?a / ?b)  \<up>  ?n = ?a  \<up>  ?n / ?b  \<up>  ?n"*)

         program = Rule.Empty_Prog}),

   		Rule.Rls_ make_rat_poly_with_parentheses,

   		Rule.Rls_ cancel_p_rls,(*FIXME:cancel_p does NOT order sometimes*)

   		Rule.Rls_ rat_reduce_1],

     program = Rule.Empty_Prog};

 (*.for simplify_Integral adapted from 'norm_Rational'.*)

 val norm_Rational_noadd_fractions = 

   Rule_Set.Sequence {id = "norm_Rational_noadd_fractions", preconds = [], 

     rew_ord = ("dummy_ord",Rewrite_Ord.function_empty), 

     asm_rls = norm_rat_erls, prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

     rules = [Rule.Rls_ discard_minus,

   		Rule.Rls_ rat_mult_poly,(* removes double fractions like a/b/c    *)

   		Rule.Rls_ make_rat_poly_with_parentheses, (*WN0510 also in(#)below*)

   		Rule.Rls_ cancel_p_rls, (*FIXME.MG:cancel_p does NOT order sometim*)

   		Rule.Rls_ norm_Rational_rls_noadd_fractions,(* the main rls (#)   *)

   		Rule.Rls_ discard_parentheses1], (* mult only                       *)

     program = Rule.Empty_Prog};

 (*.simplify terms before and after Integration such that  

    ..a.x^2/2 + b.x^3/3.. is made to ..a/2.x^2 + b/3.x^3.. (and NO

    common denominator as done by norm_Rational or make_ratpoly_in.

    This is a copy from 'make_ratpoly_in' with respective reduction of rules and

    *1* expand the term, ie. distribute * and / over +

 .*)

 val separate_bdv2 =

    Rule_Set.append_rules "separate_bdv2" collect_bdv [

     \<^rule_thm>\<open>separate_bdv\<close>, (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*)

 		\<^rule_thm>\<open>separate_bdv_n\<close>,

 		\<^rule_thm>\<open>separate_1_bdv\<close>, (*"?bdv / ?b = (1 / ?b) * ?bdv"*)

 		\<^rule_thm>\<open>separate_1_bdv_n\<close> (*"?bdv  \<up>  ?n / ?b = 1 / ?b * ?bdv  \<up>  ?n"*)

     (*

 		rule_thm>\<open>add_divide_distrib\<close> (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)*)

 		];

 val simplify_Integral = 

   Rule_Set.Sequence {id = "simplify_Integral", preconds = []:term list, 

     rew_ord = ("dummy_ord", Rewrite_Ord.function_empty),

     asm_rls = Atools_erls, prog_rls = Rule_Set.Empty,

     calc = [],  errpatts = [],

     rules = [

       \<^rule_thm>\<open>distrib_right\<close>, (*"(?z1.0 + ?z2.0) * ?w = ?z1.0 * ?w + ?z2.0 * ?w"*)

 	    \<^rule_thm>\<open>add_divide_distrib\<close>, (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)

 	     (*^^^^^ *1* ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^*)

 	    Rule.Rls_ norm_Rational_noadd_fractions,

 	    Rule.Rls_ order_add_mult_in,

 	    Rule.Rls_ discard_parentheses,

 	    (*Rule.Rls_ collect_bdv, from make_polynomial_in*)

 	    Rule.Rls_ separate_bdv2,

 	    \<^rule_eval>\<open>divide\<close> (Prog_Expr.eval_cancel "#divide_e")],

     program = Rule.Empty_Prog};      

 val integration =

   Rule_Set.Sequence {

      id="integration", preconds = [], rew_ord = ("termlessI",termlessI), 

   	 asm_rls = Rule_Def.Repeat {

        id="conditions_in_integration",  preconds = [], rew_ord = ("termlessI",termlessI), 

   		 asm_rls = Rule_Set.Empty, prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

   		 rules = [], program = Rule.Empty_Prog}, 

   	 prog_rls = Rule_Set.Empty, calc = [], errpatts = [],

   	 rules = [

       Rule.Rls_ integration_rules,

   		 Rule.Rls_ add_new_c,

   		 Rule.Rls_ simplify_Integral],

   	 program = Rule.Empty_Prog};

 val prep_rls' = Auto_Prog.prep_rls @{theory};

 \<close>

 rule_set_knowledge

   integration_rules = \<open>prep_rls' integration_rules\<close> and

   add_new_c = \<open>prep_rls' add_new_c\<close> and

   simplify_Integral = \<open>prep_rls' simplify_Integral\<close> and

   integration = \<open>prep_rls' integration\<close> and

   separate_bdv2 = \<open>prep_rls' separate_bdv2\<close> and

   norm_Rational_noadd_fractions = \<open>prep_rls' norm_Rational_noadd_fractions\<close> and

   norm_Rational_rls_noadd_fractions = \<open>prep_rls' norm_Rational_rls_noadd_fractions\<close>

 (** problems **)

 problem pbl_fun_integ : "integrate/function" =

   \<open>Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)]\<close>

   Method_Ref: "diff/integration"

   CAS: "Integrate (f_f, v_v)"

   Given: "functionTerm f_f" "integrateBy v_v"

   Find: "antiDerivative F_F"

 problem pbl_fun_integ_nam : "named/integrate/function" =

   (*here "named" is used differently from Differentiation"*)

   \<open>Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)]\<close>

   Method_Ref: "diff/integration/named"

   CAS: "Integrate (f_f, v_v)"

   Given: "functionTerm f_f" "integrateBy v_v"

   Find: "antiDerivativeName F_F"

 (** methods **)

 partial_function (tailrec) integrate :: "real \<Rightarrow> real \<Rightarrow> real"

   where

 "integrate f_f v_v = (

   let

     t_t = Take (Integral f_f D v_v)

   in

     (Rewrite_Set_Inst [(''bdv'', v_v)] ''integration'') t_t)"

 method met_diffint : "diff/integration" =

   \<open>{rew_ord="tless_true", rls'=Atools_erls, calc = [], prog_rls = Rule_Set.empty, where_rls=Rule_Set.empty,

 	  errpats = [], rew_rls = Rule_Set.empty}\<close>

   Program: integrate.simps

   Given: "functionTerm f_f" "integrateBy v_v"

   Find: "antiDerivative F_F"

 partial_function (tailrec) intergrate_named :: "real \<Rightarrow> real \<Rightarrow> (real \<Rightarrow> real) \<Rightarrow> bool"

   where

 "intergrate_named f_f v_v F_F =(

   let

     t_t = Take (F_F v_v = Integral f_f D v_v)

   in (

     (Try (Rewrite_Set_Inst [(''bdv'', v_v)] ''simplify_Integral'')) #>

     (Rewrite_Set_Inst [(''bdv'', v_v)] ''integration'')

     ) t_t)"

 method met_diffint_named : "diff/integration/named" =

   \<open>{rew_ord="tless_true", rls'=Atools_erls, calc = [], prog_rls = Rule_Set.empty, where_rls=Rule_Set.empty,

     errpats = [], rew_rls = Rule_Set.empty}\<close>

   Program: intergrate_named.simps

   Given: "functionTerm f_f" "integrateBy v_v"

   Find: "antiDerivativeName F_F"

 ML \<open>

 \<close> ML \<open>

 \<close>

 end