Theory Incomplete_Gamma.Safe_Power

section ‹A ``safe'' power operator for real and complex numbers›
theory Safe_Power
  imports "HOL-Complex_Analysis.Complex_Analysis"
begin

text ‹
  Isabelle/HOL currently has three different power operators:

     constpower, denoted termx ^ n, where n› is a natural number and x› is an element of a
      multiplicatively written monoid (classmonoid_mult).

     constpower_int, denoted termx powi n, where n› is an integer and x› is an element of
      a multiplicatively written monoid with some kind of inverse operation (typically a 
      classdivision_ring).

     constpowr, denoted termx powr y, where x› and y› are complete normed real algebra
      with a 1 and some kind of natural logarithm. In practice, the only reasonable examples of
      such a structure are probably typreal and typcomplex. It is defined as 
      propx powr y = exp (ln x * y), except when x = 0›, in which case it returns 0› by
      convention.

   All three of these operators are required, since none of them is a generalisation of another.
   For example, the constpower operator is the most restrictive in what the right argument can
   be, but it is the most general in terms of what the left argument can be.

   For the most part, the three operators agree in all the cases where they are simultaneously
   defined, but there are some caveats.

   First of all, note that different conventions apply for $0^0$: We have propx ^ n = 1 and
   propx powr n = 1 for any x›, including for x = 0›, whereas prop0 powr y = 0 for any y›, 
   including y = 0›.

   Second, for negative real x›, the constpowr operator inherits somewhat strange behaviour from
   the real-valued constln operator, which in Isabelle/HOL is defined symmetrically,
   i.e.\ propln (-x) = ln x, so we also have prop(-x) powr y = x powr y for real x› and y›. 
   This means that while prop(-1::real) ^ 3 = -1 as expected, we have 
   prop(-1::real) powr 3 = 1. It is therefore better to consider termx powr y to be undefined
   for negative real x›.

   In the following, we will define a ``safe'' version of the constpowr operator that coincides
   with constpower and constpower_int whenever applicable.
›


text ‹
  First, we define the inverse of the canonical ring homomorphism $\mathbb{Z}\to R$ using
  the choice operator, i.e.\ for any $x\in\mathbb{Z}$, the operator the_int› gives us an
  integer $n$ such that $x = n = \underbrace{1 + \ldots + 1}_{n\ \text{times}}$. 
  For $x\notin\mathbb{Z}$, the operator's behaviour is unspecified.
›
definition the_int :: "'a :: ring_char_0  int" where
  "the_int x = (THE n. x = of_int n)"

lemma the_int_of_int [simp]: "the_int (of_int n :: 'a :: ring_char_0) = n"
  unfolding the_int_def using theI'[of "λm. of_int n = (of_int m :: 'a)"] by simp

lemma the_int_eqI: "of_int n = x  the_int x = n"
  by (metis the_int_of_int)

lemma the_int_0 [simp]: "the_int 0 = 0"
  using the_int_of_int[of 0] by (simp del: the_int_of_int)

lemma the_int_1 [simp]: "the_int 1 = 1"
  using the_int_of_int[of 1] by (simp del: the_int_of_int)

lemma the_int_of_nat [simp]: "the_int (of_nat n) = int n"
  using the_int_of_int[of "int n"] by (simp del: the_int_of_int)

lemma the_int_numeral [simp]: "the_int (numeral n) = numeral n"
  using the_int_of_int[of "numeral n"] by (simp del: the_int_of_int)

lemma the_int_uminus [simp]: "x    the_int (-x) = -the_int x"
  by (metis Ints_cases of_int_minus the_int_of_int)

lemma of_int_the_int: "x    of_int (the_int x) = x"
  by (elim Ints_cases) auto

lemma the_int_add: "x    y    the_int (x + y) = the_int x + the_int y"
  by (rule the_int_eqI) (auto simp: of_int_the_int)

lemma the_int_diff: "x    y    the_int (x - y) = the_int x - the_int y"
  by (rule the_int_eqI) (auto simp: of_int_the_int)

lemma the_int_mult: "x    y    the_int (x * y) = the_int x * the_int y"
  by (rule the_int_eqI) (auto simp: of_int_the_int)

lemma the_int_power: "x    the_int (x ^ n) = the_int x ^ n"
  by (rule the_int_eqI) (auto simp: of_int_the_int)


text ‹
  Now we simply define the powr'› operator as a version of the constpowr operator that
  falls back to the constpower_int operator whenever possible.
›
definition powr' :: "'a :: {division_ring, ln}  'a  'a" (infixr powr'' 80) where
  "x powr' y = (if y   then x powi (the_int y) else x powr y)"

lemma powr'_powr: "y    x powr' y = x powr y"
  by (simp add: powr'_def)

lemma powr'_0_left [simp]: "y  0  0 powr' y = 0"
  by (auto simp: powr'_def elim!: Ints_cases)

lemma powr'_0_left_if: "0 powr' y = (if y = 0 then 1 else 0)"
  by (auto simp: powr'_def elim!: Ints_cases)

lemma powr'_of_int [simp]: "x powr' of_int n = x powi n"
  by (simp add: powr'_def)

lemma powr'_of_nat [simp]: "x powr' of_nat n = x ^ n"
  by (simp add: powr'_def)

lemma powr'_0 [simp]: "x powr' 0 = 1"
  by (simp add: powr'_def)

lemma powr'_1 [simp]: "x powr' 1 = x"
  by (simp add: powr'_def)

lemma powr'_numeral [simp]: "x powr' numeral n = x ^ numeral n"
  by (simp add: powr'_def)

text ‹
  If $x$ is positive or if $x$ is non-negative and $y$ is positive, the safe power operator
  and the normal power operator agree..
›
lemma powr'_real: "x  0  x  0  y > 0  x powr' y = x powr (y :: real)"
  by (auto simp: powr'_def powr_real_of_int' elim!: Ints_cases)

lemma powr'_real_pos: "x > 0  x powr' y = x powr (y :: real)"
  by (auto simp: powr'_def powr_real_of_int' elim!: Ints_cases)

text ‹
  For complex inputs, the two operators always agree except in the case of $0^0$.
›
lemma powr'_complex: "x  0  y  0  x powr' y = x powr (y :: complex)"
  by (auto simp: powr'_def complex_powr_of_int elim!: Ints_cases)

lemma powr'_complex_of_real:
  "x  0  y    complex_of_real x powr' of_real y = (of_real (x powr' y))"
  by (auto simp: powr'_def powr_Reals_eq elim!: Ints_cases)

lemma powr'_Reals_eq: "x  ; y  ; Re x  0  x powr' y = of_real (Re x powr' Re y)"
  by (cases "x = 0"; cases "y = 0")
     (auto elim!: Reals_cases simp: powr'_complex powr'_real powr_Reals_eq)

text ‹
  For this reason, the constpowr' operator is holomorphic in both inputs except for a branch
  cut along the non-positive reals.
›
lemma holomorphic_on_powr' [holomorphic_intros]:
  assumes "f holomorphic_on A" "g holomorphic_on A" "x. x  A  f x  0"
  shows   "(λx. f x powr' g x) holomorphic_on A"
proof -
  have "(λx. f x powr g x) holomorphic_on A"
    by (intro holomorphic_intros assms)
  also have "(λx. f x powr g x) holomorphic_on A  (λx. f x powr' g x) holomorphic_on A"
    by (intro holomorphic_cong refl) (subst powr'_complex, auto dest: assms(3))
  finally show ?thesis .
qed

lemma analytic_on_powr' [analytic_intros]:
  assumes "f analytic_on A" "g analytic_on A" "x. x  A  f x  0"
  shows   "(λx. f x powr' g x) analytic_on A"
proof -
  from assms(1) obtain B where B: "open B" "A  B" "f holomorphic_on B"
    using analytic_on_holomorphic by blast
  from assms(2) obtain C where C: "open C" "A  C" "g holomorphic_on C"
    using analytic_on_holomorphic by blast
  note [holomorphic_intros] = holomorphic_on_subset[OF B(3)] holomorphic_on_subset[OF C(3)]
  have "(λx. f x powr' g x) holomorphic_on ((B  C)  f -` (-0))"
    by (intro holomorphic_intros) auto
  moreover have "open ((B  C)  f -` (-0))" using B(1) C(1)
    by (intro continuous_open_preimage holomorphic_on_imp_continuous_on holomorphic_intros) auto
  moreover have "A  (B  C)  f -` (-0)"
    using assms(3) B C by auto
  ultimately show ?thesis
    using analytic_on_holomorphic by blast
qed

(* TODO: sort out version in library *)
lemma has_field_derivative_powr_complex [derivative_intros]:
  assumes "(f has_field_derivative f') (at x within A)"
  assumes "(g has_field_derivative g') (at x within A)"
  assumes "f x  (0 :: complex set)"
  shows   "((λx. f x powr g x) has_field_derivative 
             (f x powr g x * (f' / f x * g x + g' * ln (f x)))) (at x within A)"
proof -
  have "((λx. exp (ln (f x) * g x)) has_field_derivative 
          (exp (ln (f x) * g x) * (f' / f x * g x + g' * ln (f x)))) (at x within A)"
    by (auto intro!: derivative_eq_intros assms simp: field_simps)
  also have "(exp (ln (f x) * g x) * (f' / f x * g x + g' * ln (f x))) = 
             (f x powr g x * (f' / f x * g x + g' * ln (f x)))"
    by (use assms(3) in auto simp: powr_def mult_ac)
  also have "((λx. exp (ln (f x) * g x)) has_field_derivative ) (at x within A)  ?thesis"
  proof (intro has_field_derivative_cong_eventually)
    have "eventually (λz. z  -{0}) (nhds (f x))"
      by (intro eventually_nhds_in_open) (use assms(3) in auto)
    moreover have "continuous (at x within A) f"
      using assms(1) DERIV_continuous by blast
    hence "filterlim f (nhds (f x)) (at x within A)"
      by (simp add: continuous_within)
    ultimately have "eventually (λx. f x  -{0}) (at x within A)"
      by (rule eventually_compose_filterlim)
    thus "F x in at x within A. exp (ln (f x) * g x) = f x powr g x"
      by eventually_elim (auto simp: powr_def mult_ac)
  next
    show "exp (Ln (f x) * g x) = f x powr g x"
      by (use assms(3) in auto simp: powr_def mult_ac)
  qed
  finally show ?thesis .
qed

(* TODO: real *)
lemma has_field_derivative_powr'_complex [derivative_intros]:
  assumes "(f has_field_derivative f') (at x within A)"
  assumes "(g has_field_derivative g') (at x within A)"
  assumes "f x  (0 :: complex set)"
  shows   "((λx. f x powr' g x) has_field_derivative 
             (f x powr' g x * (f' / f x * g x + g' * ln (f x)))) (at x within A)"
proof -
  have "((λx. exp (ln (f x) * g x)) has_field_derivative 
          (exp (ln (f x) * g x) * (f' / f x * g x + g' * ln (f x)))) (at x within A)"
    by (auto intro!: derivative_eq_intros assms simp: field_simps)
  also have "(exp (ln (f x) * g x) * (f' / f x * g x + g' * ln (f x))) = 
             (f x powr' g x * (f' / f x * g x + g' * ln (f x)))"
    by (subst powr'_complex) (use assms(3) in auto simp: powr_def mult_ac)
  also have "((λx. exp (ln (f x) * g x)) has_field_derivative ) (at x within A)  ?thesis"
  proof (intro has_field_derivative_cong_eventually)
    have "eventually (λz. z  -{0}) (nhds (f x))"
      by (intro eventually_nhds_in_open) (use assms(3) in auto)
    moreover have "continuous (at x within A) f"
      using assms(1) DERIV_continuous by blast
    hence "filterlim f (nhds (f x)) (at x within A)"
      by (simp add: continuous_within)
    ultimately have "eventually (λx. f x  -{0}) (at x within A)"
      by (rule eventually_compose_filterlim)
    thus "F x in at x within A. exp (ln (f x) * g x) = f x powr' g x"
      by eventually_elim (auto simp: powr'_complex powr_def mult_ac)
  next
    show "exp (Ln (f x) * g x) = f x powr' g x"
      by (subst powr'_complex) (use assms(3) in auto simp: powr_def mult_ac)
  qed
  finally show ?thesis .
qed

lemma has_field_derivative_powr'_Ints:
  assumes "(f has_field_derivative f') (at x within A)"
  assumes "c  ( :: 'a :: {real_normed_field, ln} set)" "c    f x  0"
  shows   "((λx. f x powr' c) has_field_derivative (c * f x powr' (c-1) * f')) (at x within A)"
proof -
  from assms(2) obtain n where c: "c = of_int n"
    by (elim Ints_cases)
  have "n  0  f x  0"
    using assms(3) c
    by (metis Nats_cases dual_order.refl int_nat_eq nat_int of_int_eq_iff of_int_of_nat_eq)
  hence "((λx. f x powi n) has_field_derivative (of_int n * f x powi (n - 1) * f')) (at x within A)"
    by (auto intro!: derivative_eq_intros assms(1))
  also have "of_int n * f x powi (n - 1) * f' = c * f x powr' (of_int (n - 1)) * f'"
    by (subst powr'_of_int) (use c in auto)
  also have "of_int (n - 1) = c - 1"
    using c by simp
  finally show ?thesis
    unfolding c by (subst powr'_of_int) auto
qed

lemma continuous_on_powr'_complex [continuous_intros]:
  assumes "A  {z. Re (f z)  0  Im (f z)  0}"
  assumes "z. z  A  f z = 0  Re (g z) > 0"
  assumes "continuous_on A f" "continuous_on A g"
  shows   "continuous_on A (λz. f z powr' g z)"
proof -
  have "continuous_on A (λz. f z powr g z)"
    by (intro continuous_intros assms)
  also have "?this  ?thesis"
  proof (rule continuous_on_cong)
    fix x assume "x  A"
    hence "f x  0  g x  0"
      using assms(2)[of x] assms(1) by (auto simp: complex_eq_iff)
    thus "f x powr g x = f x powr' g x"
      by (simp add: powr'_complex)
  qed auto
  finally show ?thesis .
qed

lemma tendsto_powr'_complex [tendsto_intros]:
  fixes f g :: "_  complex"
  assumes "a  0  (a = 0  Re b > 0)" and "(f  a) F" "(g  b) F"
  shows   "((λz. f z powr' g z)  a powr' b) F"
proof -
  have "((λz. f z powr g z)  a powr b) F"
    by (rule tendsto_intros) (use assms in auto)
  also have "a  0  b  0"
    using assms(1) by auto
  hence "a powr b = a powr' b"
    by (simp add: powr'_complex)
  also have "eventually (λz. f z  0  g z  0) F"
  proof (cases "a = 0")
    case True
    have "eventually (λz. g z  -{0}) F"
      by (rule eventually_compose_filterlim[OF eventually_nhds_in_open assms(3)])
         (use True assms(1) in auto)
    thus ?thesis
      by eventually_elim auto
  next
    case False
    have "eventually (λz. f z  -{0}) F"
      by (rule eventually_compose_filterlim[OF eventually_nhds_in_open assms(2)]) (use False in auto)
    thus ?thesis
      by eventually_elim auto
  qed
  hence "eventually (λz. f z powr g z = f z powr' g z) F"
    by eventually_elim (auto simp: powr'_complex)
  hence "((λz. f z powr g z)  a powr' b) F  ?thesis"
    by (intro filterlim_cong) auto
  finally show ?thesis .
qed

(* TODO: replace version in library? *)
lemma ln_at_0': "filterlim (ln :: real  real) at_bot (at 0)"
proof -
  have "filterlim abs (at_right 0) (at (0::real))"
  proof (rule tendsto_imp_filterlim_at_right)
    show "abs 0 (0::real)"
      by (rule tendsto_rabs_zero) auto
  next
    show "eventually (λx. ¦x¦ > (0::real)) (at 0)"
      using eventually_neq_at_within by eventually_elim auto
  qed
  hence "filterlim (λx. ln ¦x¦ :: real) at_bot (at 0)"
    by (rule filterlim_compose[OF ln_at_0])
  also have "?this  filterlim (λx. ln x :: real) at_bot (at 0)"
    by (simp add: ln_real_def cong: if_cong)
  finally show ?thesis .
qed

(* TODO: replace version in library *)
lemma tendsto_powr_real [tendsto_intros]:
  fixes f g :: "_  real"
  assumes "(f  a) F" "(g  b) F"
  assumes "a = 0  b > 0" 
  shows   "((λz. f z powr g z)  a powr b) F"
proof (cases "a = 0")
  case False
  thus ?thesis
    using assms by (auto intro!: tendsto_intros)
next
  case [simp]: True
  have "((λz. if f z = 0 then 0 else exp (g z * ln (f z)))  0) F"
  proof (rule filterlim_If)
    show "((λz. exp (g z * ln (f z)))  0) (inf F (principal {z. f z  0}))"
    proof (rule filterlim_compose[OF exp_at_bot])
      show "filterlim (λz. g z * ln (f z)) at_bot (inf F (principal {z. f z  0}))"
      proof (rule filterlim_tendsto_pos_mult_at_bot)
        show "(g  b) (inf F (principal {z. f z  0}))"
          using assms(2) by (rule filterlim_mono) auto
      next
        show "filterlim (λx. ln (f x)) at_bot (inf F (principal {z. f z  0}))"
        proof (rule filterlim_compose[OF ln_at_0'])
          show "filterlim f (at 0) (inf F (principal {z. f z  0}))"
          proof (rule filterlim_atI)
            have "F x in principal {z. f z  0}. f x  0"
              by (auto simp: eventually_principal)
            thus "F x in inf F (principal {z. f z  0}). f x  0"
              by (rule filter_leD[rotated]) auto
          next
            show "(f  0) (inf F (principal {z. f z  0}))"
              using assms(1) by (rule filterlim_mono) auto
          qed
        qed
      qed (use assms(3) in auto)
    qed
  qed auto
  thus ?thesis
    by (simp add: powr_def)
qed

lemmas [tendsto_intros del] = tendsto_powr tendsto_powr'

lemma tendsto_powr'_real [tendsto_intros]:
  fixes f g :: "_  real"
  assumes "(f  a) F" "(g  b) F"
  assumes "a > 0  (a = 0  b > 0)" 
  shows   "((λz. f z powr' g z)  a powr' b) F"
proof (cases "a = 0")
  case False
  with assms have "a > 0"
    by auto
  have ev: "eventually (λx. f x  {0<..}) F"
    by (rule eventually_compose_filterlim[OF eventually_nhds_in_open assms(1)]) 
       (use a > 0 in auto)
  have "((λz. f z powr g z)  a powr b) F"
    by (rule tendsto_powr_real) (use assms in auto)
  also have "?this  ((λz. f z powr' g z)  a powr' b) F"
    using a > 0 by (intro filterlim_cong) (auto simp: powr'_real intro!: eventually_mono[OF ev])
  finally show ?thesis .
next
  case [simp]: True
  with assms have "b > 0"
    by auto
  have "((λz. f z powr' g z)  0) F"
    unfolding powr'_def
  proof (rule filterlim_If)
    have "((λz. f z powr g z)  0) F"
      by (rule tendsto_eq_intros assms(1,2))+ (use b > 0 in auto)
    thus "((λz. f z powr g z)  0) (inf F (principal {z. g z  }))"
      by (rule filterlim_mono) auto
  next
    show "((λz. f z powi the_int (g z))  0) (inf F (principal {z. g z  }))"
    proof (cases "inf F (principal {z. g z  }) = bot")
      case False
      define n' where "n' = round b"
      define n where "n = nat n'"
      have "b = n'"
      proof -
        from False have "F x in F. g x  "
          by (auto simp: trivial_limit_def eventually_inf_principal frequently_def)
        moreover have "eventually (λx. g x  ball b (1/2)) F"
          by (rule eventually_compose_filterlim[OF eventually_nhds_in_open assms(2)]) auto
        hence "eventually (λx. dist (g x) n' < 1) F"
        proof eventually_elim
          case (elim x)
          have "dist (g x) n'  dist (g x) b + dist b n'"
            by (rule dist_triangle)
          also have "dist (g x) b < 1/2"
            using elim by (simp add: dist_commute)
          also have "dist b n'  1 / 2"
            using of_int_round_abs_le[of b] by (auto simp: n'_def dist_norm abs_minus_commute)
          finally show ?case
            by simp
        qed
        hence "eventually (λx. g x    g x = n') F"
          by eventually_elim (auto elim!: Ints_cases simp: dist_of_int)
        ultimately have "F x in F. g x = n'"
          by (rule frequently_rev_mp)
        moreover have "F  bot"
          using _  bot by auto
        ultimately show "b = n'"
          using assms(2) limit_frequently_eq F  bot by blast
      qed

      have "b = n" "n > 0"
        using b = n' assms by (auto simp: n_def)

      have ev: "eventually (λx. g x = n) (inf F (principal {z. g z  }))"
      proof -
        have "eventually (λx. g x  ball b 1) F"
          by (rule eventually_compose_filterlim[OF eventually_nhds_in_open assms(2)]) auto
        hence "eventually (λx. g x    g x = n') F"
          by eventually_elim (auto simp: b = n' dist_of_int elim!: Ints_cases)
        thus ?thesis
          by (simp add: eventually_inf_principal b = n' flip: b = n)
      qed

      have "((λz. f z ^ n)  0) F"
        using assms(1) by (rule tendsto_eq_intros) (use n > 0 in auto)
      hence "((λz. f z ^ n)  0) (inf F (principal {z. g z  }))"
        by (rule filterlim_mono) auto
      also have "?this  ((λz. f z powi the_int (g z))  0) (inf F (principal {z. g z  }))"
      proof (rule filterlim_cong)
        show "F x in inf F (principal {z. g z  }). f x ^ n = f x powi the_int (g x)"
          using ev by eventually_elim auto
      qed auto
      finally show ?thesis .
    qed auto
  qed
  thus ?thesis
    using b > 0 by simp
qed

lemma continuous_powr'_complex [continuous_intros]: 
  assumes "continuous F f" "continuous F g"
  assumes "Re (f (netlimit F))  0  Im (f (netlimit F))  0"
  assumes "f (netlimit F) = 0  Re (g (netlimit F)) > 0"
  shows   "continuous F (λz. f z powr' g z :: complex)"
  using assms unfolding continuous_def
  by (intro tendsto_intros) (auto simp: complex_nonpos_Reals_iff complex_eq_iff)

end