(* 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)

 (*new'_c	      :: "real => real"        ("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^^^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'            *)

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

   integral_var:      "Integral bdv D bdv = bdv ^^^ 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 ^^^ n D bdv = bdv ^^^ (n+1) / (n + 1)"

 ML \<open>

 val thy = @{theory};

 (** eval functions **)

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

 (*.create a new unique variable 'c..' in a term; for use by Rule.Num_Calc in a rls;

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

    in the script; this will be possible if currying doesnt 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 id_of) var of

 		"c"::[] => true

 	      |	"c"::"_"::is => (case (TermC.int_of_str_opt o implode) is of

 				     SOME _ => true

 				   | NONE => false)

               | _ => false;

 	fun get_coeff c = case (Symbol.explode o id_of) c of

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

 			    | _ => 0;

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

     in 

 	case cs of

 	    [] => c

 	  | [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

 (*("new_c", ("Integrate.new'_c", eval_new_c "#new_c_"))*)

 fun eval_new_c _ _ (p as (Const ("Integrate.new'_c",_) $ t)) _ =

      SOME ((Rule.term2str p) ^ " = " ^ Rule.term2str (new_c p),

 	  Trueprop $ (mk_equality (p, new_c p)))

   | eval_new_c _ _ _ _ = NONE;

 *)

 (*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 (_:theory) =

     let val p' = case p of

 		     Const ("HOL.eq", T) $ lh $ rh => 

 		     Const ("HOL.eq", T) $ lh $ TermC.mk_add rh (new_c rh)

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

     in SOME ((Rule.term2str p) ^ " = " ^ Rule.term2str 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 ("Integrate.is'_f'_x", _)

 					   $ arg)) _ =

     if TermC.is_f_x arg

     then SOME ((Rule.term2str p) ^ " = True",

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

     else SOME ((Rule.term2str p) ^ " = False",

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

   | eval_is_f_x _ _ _ _ = NONE;

 \<close>

 setup \<open>KEStore_Elems.add_calcs

   [("add_new_c", ("Integrate.add'_new'_c", eval_add_new_c "add_new_c_")),

     ("is_f_x", ("Integrate.is'_f'_x", 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), 

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

 		     preconds = [], 

 		     rew_ord = ("termlessI",termlessI), 

 		     erls = Rule_Set.Empty, 

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

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

 			      Rule.Num_Calc ("Prog_Expr.occurs'_in", Prog_Expr.eval_occurs_in "#occurs_in_"),

 			      Rule.Thm ("not_true", TermC.num_str @{thm not_true}),

 			      Rule.Thm ("not_false",@{thm not_false})

 			      ],

 		     scr = Rule.EmptyScr}, 

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

 	 rules = [

 		  Rule.Thm ("integral_const", TermC.num_str @{thm integral_const}),

 		  Rule.Thm ("integral_var", TermC.num_str @{thm integral_var}),

 		  Rule.Thm ("integral_add", TermC.num_str @{thm integral_add}),

 		  Rule.Thm ("integral_mult", TermC.num_str @{thm integral_mult}),

 		  Rule.Thm ("integral_pow", TermC.num_str @{thm integral_pow}),

 		  Rule.Num_Calc ("Groups.plus_class.plus", (**)eval_binop "#add_")(*for n+1*)

 		  ],

 	 scr = Rule.EmptyScr};

 \<close>

 ML \<open>

 val add_new_c = 

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

 	 rew_ord = ("termlessI",termlessI), 

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

 		     preconds = [], 

 		     rew_ord = ("termlessI",termlessI), 

 		     erls = Rule_Set.Empty, 

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

 		     rules = [Rule.Num_Calc ("Prog_Expr.matches", Prog_Expr.eval_matches""),

 			      Rule.Num_Calc ("Integrate.is'_f'_x", 

 				    eval_is_f_x "is_f_x_"),

 			      Rule.Thm ("not_true", TermC.num_str @{thm not_true}),

 			      Rule.Thm ("not_false", TermC.num_str @{thm not_false})

 			      ],

 		     scr = Rule.EmptyScr}, 

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

 	 rules = [ (*Rule.Thm ("call_for_new_c", TermC.num_str @{thm call_for_new_c}),*)

 		   Rule.Cal1 ("Integrate.add'_new'_c", eval_add_new_c "new_c_")

 		   ],

 	 scr = Rule.EmptyScr};

 \<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",Rule.dummy_ord), 

      erls = norm_rat_erls, srls = 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",Rule.dummy_ord), 

 		       erls = (*FIXME.WN051028 Rule_Set.empty,*)

 		       Rule_Set.append_rules "Rule_Set.empty-is_polyexp" Rule_Set.empty

 				  [Rule.Num_Calc ("Poly.is'_polyexp", 

 					 eval_is_polyexp "")],

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

 				  rules = [Rule.Thm ("rat_mult", TermC.num_str @{thm rat_mult}),

 	       (*"?a / ?b * (?c / ?d) = ?a * ?c / (?b * ?d)"*)

 	       Rule.Thm ("rat_mult_poly_l", TermC.num_str @{thm rat_mult_poly_l}),

 	       (*"?c is_polyexp ==> ?c * (?a / ?b) = ?c * ?a / ?b"*)

 	       Rule.Thm ("rat_mult_poly_r", TermC.num_str @{thm rat_mult_poly_r}),

 	       (*"?c is_polyexp ==> ?a / ?b * ?c = ?a * ?c / ?b"*)

 	       Rule.Thm ("real_divide_divide1_mg",

                      TermC.num_str @{thm real_divide_divide1_mg}),

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

 	       Rule.Thm ("divide_divide_eq_right", 

                      TermC.num_str @{thm divide_divide_eq_right}),

 	       (*"?x / (?y / ?z) = ?x * ?z / ?y"*)

 	       Rule.Thm ("divide_divide_eq_left",

                      TermC.num_str @{thm divide_divide_eq_left}),

 	       (*"?x / ?y / ?z = ?x / (?y * ?z)"*)

 	       Rule.Num_Calc ("Rings.divide_class.divide", Prog_Expr.eval_cancel "#divide_e"),

 	       Rule.Thm ("rat_power", TermC.num_str @{thm rat_power})

 		(*"(?a / ?b) ^^^ ?n = ?a ^^^ ?n / ?b ^^^ ?n"*)

 	       ],

       scr = Rule.EmptyScr

       }),

 		Rule.Rls_ make_rat_poly_with_parentheses,

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

 		Rule.Rls_ rat_reduce_1

 		],

        scr = Rule.EmptyScr

        };

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

 val norm_Rational_noadd_fractions = 

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

        rew_ord = ("dummy_ord",Rule.dummy_ord), 

        erls = norm_rat_erls, srls = 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                       *)

 		],

        scr = Rule.EmptyScr

        };

 (*.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 ("separate_bdv", TermC.num_str @{thm separate_bdv}),

 		(*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*)

 		Rule.Thm ("separate_bdv_n", TermC.num_str @{thm separate_bdv_n}),

 		Rule.Thm ("separate_1_bdv",  TermC.num_str @{thm separate_1_bdv}),

 		(*"?bdv / ?b = (1 / ?b) * ?bdv"*)

 		Rule.Thm ("separate_1_bdv_n",  TermC.num_str @{thm separate_1_bdv_n})(*,

 			  (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*)

 			  *****Rule.Thm ("add_divide_distrib", 

 			  ***** TermC.num_str @{thm add_divide_distrib})

 			  (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)----------*)

 		];

 val simplify_Integral = 

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

        rew_ord = ("dummy_ord", Rule.dummy_ord),

       erls = Atools_erls, srls = Rule_Set.Empty,

       calc = [],  errpatts = [],

       rules = [Rule.Thm ("distrib_right", TermC.num_str @{thm distrib_right}),

  	       (*"(?z1.0 + ?z2.0) * ?w = ?z1.0 * ?w + ?z2.0 * ?w"*)

 	       Rule.Thm ("add_divide_distrib", TermC.num_str @{thm add_divide_distrib}),

  	       (*"(?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.Num_Calc ("Rings.divide_class.divide", Prog_Expr.eval_cancel "#divide_e")

 	       ],

       scr = Rule.EmptyScr};      

 (*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_polynomial_in' with insertions from 

    'make_ratpoly_in' 

 THIS IS KEPT FOR COMPARISON ............................................   

 * val simplify_Integral = prep_rls'(

 *   Rule_Set.Seqence {id = "", preconds = []:term list, 

 *        rew_ord = ("dummy_ord", Rule.dummy_ord),

 *       erls = Atools_erls, srls = Rule_Set.Empty,

 *       calc = [], (*asm_thm = [],*)

 *       rules = [Rule.Rls_ expand_poly,

 * 	       Rule.Rls_ order_add_mult_in,

 * 	       Rule.Rls_ simplify_power,

 * 	       Rule.Rls_ collect_numerals,

 * 	       Rule.Rls_ reduce_012,

 * 	       Rule.Thm ("realpow_oneI", TermC.num_str @{thm realpow_oneI}),

 * 	       Rule.Rls_ discard_parentheses,

 * 	       Rule.Rls_ collect_bdv,

 * 	       (*below inserted from 'make_ratpoly_in'*)

 * 	       Rule.Rls_ (Rule_Set.append_rules "separate_bdv"

 * 			 collect_bdv

 * 			 [Rule.Thm ("separate_bdv", TermC.num_str @{thm separate_bdv}),

 * 			  (*"?a * ?bdv / ?b = ?a / ?b * ?bdv"*)

 * 			  Rule.Thm ("separate_bdv_n", TermC.num_str @{thm separate_bdv_n}),

 * 			  Rule.Thm ("separate_1_bdv", TermC.num_str @{thm separate_1_bdv}),

 * 			  (*"?bdv / ?b = (1 / ?b) * ?bdv"*)

 * 			  Rule.Thm ("separate_1_bdv_n", TermC.num_str @{thm separate_1_bdv_n})(*,

 * 			  (*"?bdv ^^^ ?n / ?b = 1 / ?b * ?bdv ^^^ ?n"*)

 * 			  Rule.Thm ("add_divide_distrib", 

 * 				  TermC.num_str @{thm add_divide_distrib})

 * 			   (*"(?x + ?y) / ?z = ?x / ?z + ?y / ?z"*)*)

 * 			  ]),

 * 	       Rule.Num_Calc ("Rings.divide_class.divide"  , eval_cancel "#divide_e")

 * 	       ],

 *       scr = Rule.EmptyScr

 *       }); 

 .......................................................................*)

 val integration = 

     Rule_Set.Seqence {id="integration", preconds = [], 

 	 rew_ord = ("termlessI",termlessI), 

 	 erls = Rule_Def.Repeat {id="conditions_in_integration", 

 		     preconds = [], 

 		     rew_ord = ("termlessI",termlessI), 

 		     erls = Rule_Set.Empty, 

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

 		     rules = [],

 		     scr = Rule.EmptyScr}, 

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

 	 rules = [ Rule.Rls_ integration_rules,

 		   Rule.Rls_ add_new_c,

 		   Rule.Rls_ simplify_Integral

 		   ],

 	 scr = Rule.EmptyScr};

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

 \<close>

 setup \<open>KEStore_Elems.add_rlss 

   [("integration_rules", (Context.theory_name @{theory}, prep_rls' integration_rules)), 

   ("add_new_c", (Context.theory_name @{theory}, prep_rls' add_new_c)), 

   ("simplify_Integral", (Context.theory_name @{theory}, prep_rls' simplify_Integral)), 

   ("integration", (Context.theory_name @{theory}, prep_rls' integration)), 

   ("separate_bdv2", (Context.theory_name @{theory}, prep_rls' separate_bdv2)),

   ("norm_Rational_noadd_fractions", (Context.theory_name @{theory},

     prep_rls' norm_Rational_noadd_fractions)), 

   ("norm_Rational_rls_noadd_fractions", (Context.theory_name @{theory},

     prep_rls' norm_Rational_rls_noadd_fractions))]\<close>

 (** problems **)

 setup \<open>KEStore_Elems.add_pbts

   [(Specify.prep_pbt thy "pbl_fun_integ" [] Celem.e_pblID

       (["integrate","function"],

         [("#Given" ,["functionTerm f_f", "integrateBy v_v"]),

           ("#Find"  ,["antiDerivative F_F"])],

         Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)], 

         SOME "Integrate (f_f, v_v)", 

         [["diff","integration"]])),

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

     (Specify.prep_pbt thy "pbl_fun_integ_nam" [] Celem.e_pblID

       (["named","integrate","function"],

         [("#Given" ,["functionTerm f_f", "integrateBy v_v"]),

           ("#Find"  ,["antiDerivativeName F_F"])],

         Rule_Set.append_rules "empty" Rule_Set.empty [(*for preds in where_*)], 

         SOME "Integrate (f_f, v_v)", 

         [["diff","integration","named"]]))]\<close>

 (** 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)"

 setup \<open>KEStore_Elems.add_mets

     [Specify.prep_met thy "met_diffint" [] Celem.e_metID

 	    (["diff","integration"],

 	      [("#Given" ,["functionTerm f_f", "integrateBy v_v"]), ("#Find"  ,["antiDerivative F_F"])],

 	      {rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = Rule_Set.empty, prls=Rule_Set.empty,

 	        crls = Atools_erls, errpats = [], nrls = Rule_Set.empty},

 	      @{thm integrate.simps})]

 \<close>

 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)"

 setup \<open>KEStore_Elems.add_mets

     [Specify.prep_met thy "met_diffint_named" [] Celem.e_metID

 	    (["diff","integration","named"],

 	      [("#Given" ,["functionTerm f_f", "integrateBy v_v"]),

 	        ("#Find"  ,["antiDerivativeName F_F"])],

 	      {rew_ord'="tless_true", rls'=Atools_erls, calc = [], srls = Rule_Set.empty, prls=Rule_Set.empty,

           crls = Atools_erls, errpats = [], nrls = Rule_Set.empty},

         @{thm intergrate_named.simps})]

 \<close>

 end