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

  (*Script-names*)

  IntegrationScript      :: "[real,real,  real] => real"

                  ("((Script IntegrationScript (_ _ =))// (_))" 9)

  NamedIntegrationScript :: "[real,real, real=>real,  bool] => bool"

                  ("((Script NamedIntegrationScript (_ _ _=))// (_))" 9)

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 {*

val thy = @{theory};

(** eval functions **)

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

(*.create a new unique variable 'c..' in a term; for use by 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 (int_of_str 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 int_of_str o implode) is

			    | _ => 0;

        val cs = filter selc (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 ((term2str p) ^ " = " ^ 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 $ mk_add rh (new_c rh)

		   | p => mk_add p (new_c p)

    in SOME ((term2str p) ^ " = " ^ term2str p',

	  Trueprop $ (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 is_f_x arg

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

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

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

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

  | eval_is_f_x _ _ _ _ = NONE;

calclist':= overwritel (!calclist', 

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

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

    ]);

*}

setup {* 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_"))] *}

ML {*

(** rulesets **)

(*.rulesets for integration.*)

val integration_rules = 

    Rls {id="integration_rules", preconds = [], 

	 rew_ord = ("termlessI",termlessI), 

	 erls = Rls {id="conditions_in_integration_rules", 

		     preconds = [], 

		     rew_ord = ("termlessI",termlessI), 

		     erls = Erls, 

		     srls = Erls, calc = [], errpatts = [],

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

			      Calc ("Atools.occurs'_in", 

				    eval_occurs_in "#occurs_in_"),

			      Thm ("not_true",num_str @{thm not_true}),

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

			      ],

		     scr = EmptyScr}, 

	 srls = Erls, calc = [], errpatts = [],

	 rules = [

		  Thm ("integral_const",num_str @{thm integral_const}),

		  Thm ("integral_var",num_str @{thm integral_var}),

		  Thm ("integral_add",num_str @{thm integral_add}),

		  Thm ("integral_mult",num_str @{thm integral_mult}),

		  Thm ("integral_pow",num_str @{thm integral_pow}),

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

		  ],

	 scr = EmptyScr};

*}

ML {*

val add_new_c = 

    Seq {id="add_new_c", preconds = [], 

	 rew_ord = ("termlessI",termlessI), 

	 erls = Rls {id="conditions_in_add_new_c", 

		     preconds = [], 

		     rew_ord = ("termlessI",termlessI), 

		     erls = Erls, 

		     srls = Erls, calc = [], errpatts = [],

		     rules = [Calc ("Tools.matches", eval_matches""),

			      Calc ("Integrate.is'_f'_x", 

				    eval_is_f_x "is_f_x_"),

			      Thm ("not_true",num_str @{thm not_true}),

			      Thm ("not_false",num_str @{thm not_false})

			      ],

		     scr = EmptyScr}, 

	 srls = Erls, calc = [], errpatts = [],

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

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

		   ],

	 scr = EmptyScr};

*}

ML {*

(*.rulesets for simplifying Integrals.*)

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

val norm_Rational_rls_noadd_fractions = 

Rls {id = "norm_Rational_rls_noadd_fractions", preconds = [], 

     rew_ord = ("dummy_ord",dummy_ord), 

     erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

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

	      Rls_ (*rat_mult_div_pow original corrected WN051028*)

		  (Rls {id = "rat_mult_div_pow", preconds = [], 

		       rew_ord = ("dummy_ord",dummy_ord), 

		       erls = (*FIXME.WN051028 e_rls,*)

		       append_rls "e_rls-is_polyexp" e_rls

				  [Calc ("Poly.is'_polyexp", 

					 eval_is_polyexp "")],

				  srls = Erls, calc = [], errpatts = [],

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

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

	       Thm ("rat_mult_poly_l",num_str @{thm rat_mult_poly_l}),

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

	       Thm ("rat_mult_poly_r",num_str @{thm rat_mult_poly_r}),

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

	       Thm ("real_divide_divide1_mg",

                     num_str @{thm real_divide_divide1_mg}),

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

	       Thm ("divide_divide_eq_right", 

                     num_str @{thm divide_divide_eq_right}),

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

	       Thm ("divide_divide_eq_left",

                     num_str @{thm divide_divide_eq_left}),

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

	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e"),

	       Thm ("rat_power", num_str @{thm rat_power})

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

	       ],

      scr = Prog ((term_of o the o (parse thy)) "empty_script")

      }),

		Rls_ make_rat_poly_with_parentheses,

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

		Rls_ rat_reduce_1

		],

       scr = Prog ((term_of o the o (parse thy)) "empty_script")

       }:rls;

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

val norm_Rational_noadd_fractions = 

   Seq {id = "norm_Rational_noadd_fractions", preconds = [], 

       rew_ord = ("dummy_ord",dummy_ord), 

       erls = norm_rat_erls, srls = Erls, calc = [], errpatts = [],

       rules = [Rls_ discard_minus,

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

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

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

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

		Rls_ discard_parentheses1 (* mult only                       *)

		],

       scr = Prog ((term_of o the o (parse thy)) "empty_script")

       }:rls;

(*.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 =

    append_rls "separate_bdv2"

	       collect_bdv

	       [Thm ("separate_bdv", num_str @{thm separate_bdv}),

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

		Thm ("separate_bdv_n", num_str @{thm separate_bdv_n}),

		Thm ("separate_1_bdv", num_str @{thm separate_1_bdv}),

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

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

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

			  *****Thm ("add_divide_distrib", 

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

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

		];

val simplify_Integral = 

  Seq {id = "simplify_Integral", preconds = []:term list, 

       rew_ord = ("dummy_ord", dummy_ord),

      erls = Atools_erls, srls = Erls,

      calc = [],  errpatts = [],

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

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

	       Thm ("add_divide_distrib",num_str @{thm add_divide_distrib}),

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

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

	       Rls_ norm_Rational_noadd_fractions,

	       Rls_ order_add_mult_in,

	       Rls_ discard_parentheses,

	       (*Rls_ collect_bdv, from make_polynomial_in*)

	       Rls_ separate_bdv2,

	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")

	       ],

      scr = EmptyScr}:rls;      

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

*   Seq {id = "", preconds = []:term list, 

*        rew_ord = ("dummy_ord", dummy_ord),

*       erls = Atools_erls, srls = Erls,

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

*       rules = [Rls_ expand_poly,

* 	       Rls_ order_add_mult_in,

* 	       Rls_ simplify_power,

* 	       Rls_ collect_numerals,

* 	       Rls_ reduce_012,

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

* 	       Rls_ discard_parentheses,

* 	       Rls_ collect_bdv,

* 	       (*below inserted from 'make_ratpoly_in'*)

* 	       Rls_ (append_rls "separate_bdv"

* 			 collect_bdv

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

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

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

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

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

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

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

* 			  Thm ("add_divide_distrib", 

* 				 num_str @{thm add_divide_distrib})

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

* 			  ]),

* 	       Calc ("Fields.inverse_class.divide"  ,eval_cancel "#divide_e")

* 	       ],

*       scr = EmptyScr

*       }:rls); 

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

val integration = 

    Seq {id="integration", preconds = [], 

	 rew_ord = ("termlessI",termlessI), 

	 erls = Rls {id="conditions_in_integration", 

		     preconds = [], 

		     rew_ord = ("termlessI",termlessI), 

		     erls = Erls, 

		     srls = Erls, calc = [], errpatts = [],

		     rules = [],

		     scr = EmptyScr}, 

	 srls = Erls, calc = [], errpatts = [],

	 rules = [ Rls_ integration_rules,

		   Rls_ add_new_c,

		   Rls_ simplify_Integral

		   ],

	 scr = EmptyScr};

ruleset' := 

overwritelthy @{theory} (!ruleset', 

	    [("integration_rules", prep_rls integration_rules),

	     ("add_new_c", prep_rls add_new_c),

	     ("simplify_Integral", prep_rls simplify_Integral),

	     ("integration", prep_rls integration),

	     ("separate_bdv2", separate_bdv2),

	     ("norm_Rational_noadd_fractions", norm_Rational_noadd_fractions),

	     ("norm_Rational_rls_noadd_fractions", 

	      norm_Rational_rls_noadd_fractions)

	     ]);

*}

setup {* 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))] *}

ML {*

(** problems **)

store_pbt

 (prep_pbt thy "pbl_fun_integ" [] e_pblID

 (["integrate","function"],

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

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

  ],

  append_rls "e_rls" e_rls [(*for preds in where_*)], 

  SOME "Integrate (f_f, v_v)", 

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

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

store_pbt

 (prep_pbt thy "pbl_fun_integ_nam" [] e_pblID

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

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

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

  ],

  append_rls "e_rls" e_rls [(*for preds in where_*)], 

  SOME "Integrate (f_f, v_v)", 

  [["diff","integration","named"]]));

*}

ML {*

(** methods **)

store_met

    (prep_met thy "met_diffint" [] e_metID

	      (["diff","integration"],

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

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

		],

	       {rew_ord'="tless_true", rls'=Atools_erls, calc = [], 

		srls = e_rls, 

		prls=e_rls,

	     crls = Atools_erls, errpats = [], nrls = e_rls},

"Script IntegrationScript (f_f::real) (v_v::real) =                " ^

"  (let t_t = Take (Integral f_f D v_v)                             " ^

"   in (Rewrite_Set_Inst [(bdv,v_v)] integration False) (t_t::real))"

));

  *}

ML {*

store_met

    (prep_met thy "met_diffint_named" [] 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 = e_rls, 

		prls=e_rls,

		crls = Atools_erls, errpats = [], nrls = e_rls},

"Script NamedIntegrationScript (f_f::real) (v_v::real) (F_F::real=>real) = " ^

"  (let t_t = Take (F_F v_v = Integral f_f D v_v)                            " ^

"   in ((Try (Rewrite_Set_Inst [(bdv,v_v)] simplify_Integral False)) @@  " ^

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

    ));

*}

end