Theory Positive_Operators

section Positive_Operators› -- Positive bounded operators›

theory Positive_Operators
  imports
    Ordinary_Differential_Equations.Cones
    Misc_Tensor_Product_BO
    Strong_Operator_Topology
begin

no_notation Infinite_Set_Sum.abs_summable_on (infix "abs'_summable'_on" 50)
hide_const (open) Infinite_Set_Sum.abs_summable_on
hide_fact (open) Infinite_Set_Sum.abs_summable_on_Sigma_iff

unbundle cblinfun_notation

lemma cinner_pos_if_pos: f C (A *V f)  0 if A  0
  using less_eq_cblinfun_def that by force

definition sqrt_op :: ('a::chilbert_space CL 'a)  ('a CL 'a) where
  sqrt_op a = (if (b :: 'a CL 'a. b  0  b* oCL b = a) then (SOME b. b  0  b* oCL b = a) else 0)

lemma sqrt_op_nonpos: sqrt_op a = 0 if ¬ a  0
proof -
  have ¬ (b. b  0  b* oCL b = a)
    using positive_cblinfun_squareI that by blast
  then show ?thesis
    by (auto simp add: sqrt_op_def)
qed

lemma generalized_Cauchy_Schwarz:
  fixes inner A
  assumes Apos: A  0
  defines "inner x y  x C (A *V y)"
  shows complex_of_real ((norm (inner x y))2)  inner x x * inner y y
proof (cases inner y y = 0)
  case True
  have [simp]: inner (s *C x) y = cnj s * inner x y for s x y
    by (simp add: assms(2))
  have [simp]: inner x (s *C y) = s * inner x y for s x y
    by (simp add: assms(2) cblinfun.scaleC_right)
  have [simp]: inner (x - x') y = inner x y - inner x' y for x x' y
    by (simp add: cinner_diff_left inner_def)
  have [simp]: inner x (y - y') = inner x y - inner x y' for x y y'
    by (simp add: cblinfun.diff_right cinner_diff_right inner_def)
  have Re0: Re (inner x y) = 0 for x
  proof -
    have *: Re (inner x y) = (inner x y + inner y x) / 2
      by (smt (verit, del_insts) assms(1) assms(2) cinner_adj_left cinner_commute complex_Re_numeral complex_add_cnj field_sum_of_halves numeral_One numeral_plus_numeral of_real_divide of_real_numeral one_complex.simps(1) positive_hermitianI semiring_norm(2))
    have 0  Re (inner (x - s *C y) (x - s *C y)) for s
      by (metis Re_mono assms(1) assms(2) cinner_pos_if_pos zero_complex.simps(1))
    also have  s = Re (inner x x) - s * 2 * Re (inner x y) for s
      apply (auto simp: True)
      by (smt (verit, ccfv_threshold) Re_complex_of_real assms(1) assms(2) cinner_adj_right cinner_commute complex_add_cnj diff_minus_eq_add minus_complex.simps(1) positive_hermitianI uminus_complex.sel(1))
    finally show Re (inner x y) = 0
      by (metis add_le_same_cancel1 ge_iff_diff_ge_0 nonzero_eq_divide_eq not_numeral_le_zero zero_neq_numeral)
  qed
  have Im (inner x y) = Re (inner (imaginary_unit *C x) y)
    by simp
  also have  = 0
    by (rule Re0)
  finally have inner x y = 0
    using Re0[of x]
    using complex_eq_iff zero_complex.simps(1) zero_complex.simps(2) by presburger
  then show ?thesis
    by (auto simp: True)
next
  case False
  have inner_commute: inner x y = cnj (inner y x)
    by (metis Apos cinner_adj_left cinner_commute' inner_def positive_hermitianI)
  have [simp]: "cnj (inner y y) = inner y y" for y
    by (metis assms(1) cinner_adj_right cinner_commute' inner_def positive_hermitianI)
  define r where "r = cnj (inner x y) / inner y y"
  have "0  inner (x - scaleC r y) (x - scaleC r y)"
    by (simp add: Apos inner_def cinner_pos_if_pos)
  also have " = inner x x - r * inner x y - cnj r * inner y x + r * cnj r * inner y y"
    unfolding cinner_diff_left cinner_diff_right cinner_scaleC_left cinner_scaleC_right inner_def
    by (smt (verit, ccfv_threshold) cblinfun.diff_right cblinfun.scaleC_right cblinfun_cinner_right.rep_eq cinner_scaleC_left cinner_scaleC_right diff_add_eq diff_diff_eq2 mult.assoc)
  also have " = inner x x - inner y x * cnj r"
    unfolding r_def by auto
  also have " = inner x x - inner x y * cnj (inner x y) / inner y y"
    unfolding r_def
    by (metis assms(1) assms(2) cinner_adj_right cinner_commute complex_cnj_divide mult.commute positive_hermitianI times_divide_eq_left)
  finally have "0  inner x x - inner x y * cnj (inner x y) / inner y y" .
  hence "inner x y * cnj (inner x y) / inner y y  inner x x"
    by (simp add: le_diff_eq)
  hence (norm (inner x y)) ^ 2 / inner y y  inner x x
    using complex_norm_square by presburger
  then show ?thesis
    by (metis False assms(1) assms(2) cinner_pos_if_pos mult_right_mono nonzero_eq_divide_eq)
qed

lemma sandwich_pos[intro]: sandwich b a  0 if a  0
  by (metis (no_types, opaque_lifting) positive_cblinfunI cblinfun_apply_cblinfun_compose cinner_adj_left cinner_pos_if_pos sandwich_apply that)

lemma cblinfun_power_pos: cblinfun_power a n  0 if a  0
proof (cases even n)
  case True
  have 0  (cblinfun_power a (n div 2))* oCL (cblinfun_power a (n div 2))
    using positive_cblinfun_squareI by blast
  also have  = cblinfun_power a (n div 2 + n div 2)
    by (metis cblinfun_power_adj cblinfun_power_compose positive_hermitianI that)
  also from True have  = cblinfun_power a n
    by (metis add_self_div_2 div_plus_div_distrib_dvd_right) 
  finally show ?thesis
    by -
next
  case False
  have 0  sandwich (cblinfun_power a (n div 2)) a
    using a  0 by (rule sandwich_pos)
  also have  = cblinfun_power a (n div 2 + 1 + n div 2)
    unfolding sandwich_apply
    by (metis (no_types, lifting) One_nat_def cblinfun_compose_id_right cblinfun_power_0 cblinfun_power_Suc' cblinfun_power_adj cblinfun_power_compose positive_hermitianI that)
  also from False have  = cblinfun_power a n
    by (smt (verit, del_insts) Suc_1 add.commute add.left_commute add_mult_distrib2 add_self_div_2 nat.simps(3) nonzero_mult_div_cancel_left odd_two_times_div_two_succ)
  finally show ?thesis
    by -
qed



(* Proof follows https://link.springer.com/article/10.1007%2FBF01448052,
      @{cite wecken35linearer} *)
lemma sqrt_op_existence:
  fixes A :: 'a::chilbert_space CL 'a::chilbert_space
  assumes Apos: A  0
  shows B. B  0  B oCL B = A  (F. A oCL F = F oCL A  B oCL F = F oCL B)
           B  closure (cspan (range (cblinfun_power A)))
proof -
  define k S where k = norm A and S = A /R k - id_cblinfun
  have S  0
  proof (rule cblinfun_leI)
    fix x :: 'a assume [simp]: norm x = 1
    with assms have aux1: complex_of_real (inverse (norm A)) * (x C (A *V x))  1
      by (smt (verit, del_insts) Reals_cnj_iff cinner_adj_left cinner_commute cinner_scaleR_left cinner_scaleR_right cmod_Re complex_inner_class.Cauchy_Schwarz_ineq2 left_inverse less_eq_complex_def linordered_field_class.inverse_nonnegative_iff_nonnegative mult_cancel_left2 mult_left_mono norm_cblinfun norm_ge_zero norm_mult norm_of_real norm_one positive_hermitianI reals_zero_comparable zero_less_one_class.zero_le_one)
    show x C (S *V x)  x C (0 *V x)
      by (auto simp: S_def cinner_diff_right cblinfun.diff_left scaleR_scaleC cdot_square_norm k_def complex_of_real_mono_iff[where y=1, simplified]
          simp flip: assms of_real_inverse of_real_power of_real_mult power_mult_distrib power_inverse
          intro!: power_le_one aux1)
  qed
  have [simp]: S* = S
    using S  0 adj_0 comparable_hermitean' selfadjoint_def by blast
  have - id_cblinfun  S
    by (simp add: S_def assms k_def scaleR_nonneg_nonneg)
  then have norm (S *V f)  norm f for f
  proof -
    have 1: - S  0
      by (simp add: S  0)
    have 2: f C (- S *V f)  f C f
      by (metis - id_cblinfun  S id_cblinfun_apply less_eq_cblinfun_def minus_le_iff)

    have (norm (S *V f))^4 = complex_of_real ((cmod ((- S *V f) C (- S *V f)))2)
      apply (auto simp: power4_eq_xxxx cblinfun.minus_left complex_of_real_cmod power2_eq_square simp flip: power2_norm_eq_cinner)
      by (smt (verit, ccfv_SIG) complex_of_real_cmod mult.assoc norm_ge_zero norm_mult norm_of_real of_real_mult)
    also have   (- S *V f) C (- S *V - S *V f) * (f C (- S *V f))
      apply (rule generalized_Cauchy_Schwarz[where A=-S and x = -S *V f and y = f])
      by (fact 1)
    also have   (- S *V f) C (- S *V - S *V f) * (f C f)
      using 2 apply (rule mult_left_mono)
      using "1" cinner_pos_if_pos by blast
    also have   (- S *V f) C (- S *V f) * (f C f)
      apply (rule mult_right_mono)
      apply (metis - id_cblinfun  S id_cblinfun_apply less_eq_cblinfun_def neg_le_iff_le verit_minus_simplify(4))
      by simp
    also have  = (norm (-S *V f))2 * (norm f)2
      by (simp add: cdot_square_norm)
    also have  = (norm (S *V f))2 * (norm f)2
      by (simp add: cblinfun.minus_left)
    finally have norm (S *V f) ^ 4  (norm (S *V f))2 * (norm f)2
      using complex_of_real_mono_iff by blast
    then have (norm (S *V f))2  (norm f)2
      by (smt (verit, best) complex_of_real (norm (S *V f) ^ 4) = complex_of_real ((cmod ((- S *V f) C (- S *V f)))2) cblinfun.real.minus_left cinner_ge_zero cmod_Re mult_cancel_left mult_left_mono norm_minus_cancel of_real_eq_iff power2_eq_square power2_norm_eq_cinner' zero_less_norm_iff)
    then show norm (S *V f)  norm f
      by auto
  qed
  then have norm_Snf: norm (cblinfun_power S n *V f)  norm f for f n
    by (induction n, auto simp: cblinfun_power_Suc' intro: order.trans)
  have fSnf: cmod (f C (cblinfun_power S n *V f))  cmod (f C f) for f n
    by (smt (z3) One_nat_def Re_complex_of_real Suc_1 cdot_square_norm cinner_ge_zero cmod_Re complex_inner_class.Cauchy_Schwarz_ineq2 mult.commute mult_cancel_right1 mult_left_mono norm_Snf norm_ge_zero power_0 power_Suc)
  from norm_Snf have norm_Sn: norm (cblinfun_power S n)  1 for n
    apply (rule_tac norm_cblinfun_bound)
    by auto
  define b where b = (λn. (1/2 gchoose n) *R cblinfun_power S n)
  define B0 B where B0 = infsum b UNIV and B = sqrt k *R B0

  have sum_norm_b: (nF. norm (b n))  3 (is ?lhs  ?rhs) if finite F for F
  proof -
    have [simp]: 1 / 2 :: real = 1
      by (simp add: ceiling_eq_iff)
    from finite F obtain d where F  {..d} and [simp]: d > 0
      by (metis Icc_subset_Iic_iff atLeast0AtMost bot_nat_0.extremum bot_nat_0.not_eq_extremum dual_order.trans finite_nat_iff_bounded_le less_one)

    have ?lhs = (nF. norm ((1 / 2 gchoose n) *R (cblinfun_power S n)))
      by (simp add: b_def scaleR_cblinfun.rep_eq)
    also have   (nF. abs ((1 / 2 gchoose n)))
      apply (auto intro!: sum_mono)
      using norm_Sn
      by (metis norm_cmul_rule_thm norm_scaleR verit_prod_simplify(2))
    also have   (nd. abs (1/2 gchoose n))
      using F  {..d} by (auto intro!: mult_right_mono sum_mono2)
    also have  = (2 - (- 1) ^ d * (- (1 / 2) gchoose d))
      apply (subst gbinomial_sum_lower_abs)
      by auto
    also have   (2 + norm (- (1/2) gchoose d :: real))
      apply (auto intro!: mult_right_mono)
      by (smt (verit) left_minus_one_mult_self mult.assoc mult_minus_left power2_eq_iff power2_eq_square)
    also have   3
      apply (subgoal_tac abs (- (1/2) gchoose d :: real)  1)
      apply (metis add_le_cancel_left is_num_normalize(1) mult.commute mult_left_mono norm_ge_zero numeral_Bit0 numeral_Bit1 one_add_one real_norm_def)
      apply (rule abs_gbinomial_leq1)
      by auto
    finally show ?thesis
      by -
  qed

  have has_sum_b: (b has_sum B0) UNIV
    apply (auto intro!: has_sum_infsum abs_summable_summable[where f=b] bdd_aboveI[where M=3] simp: B0_def abs_summable_iff_bdd_above)
    using sum_norm_b
    by simp

  have B0  0
  proof (rule positive_cblinfunI)
    fix f :: 'a assume [simp]: norm f = 1
    from has_sum_b
    have sum1: (λn. f C (b n *V f)) summable_on UNIV
      apply (intro summable_on_cinner_left summable_on_cblinfun_apply_left)
      by (simp add: has_sum_imp_summable)
    have sum2: (λx. - (complex_of_real ¦1 / 2 gchoose x¦ * (f C f))) summable_on UNIV - {0}
      apply (rule abs_summable_summable)
      using gbinomial_abs_summable_1[of 1/2]
      by (auto simp add: cnorm_eq_1[THEN iffD1])
    from sum1 have sum3: (λn. complex_of_real (1 / 2 gchoose n) * (f C (cblinfun_power S n *V f))) summable_on UNIV - {0}
      unfolding b_def
      by (metis (no_types, lifting) cinner_scaleR_right finite.emptyI finite_insert 
          scaleR_cblinfun.rep_eq summable_on_cofin_subset summable_on_cong)

    have aux: a  - b if norm a  norm b and a   and b  0 for a b :: complex
      using cmod_eq_Re complex_is_Real_iff less_eq_complex_def that(1) that(2) that(3) by force

    from has_sum_b
    have f C (B0 *V f) = (n. f C (b n *V f))
      by (metis B0_def infsum_cblinfun_apply_left infsum_cinner_left summable_on_cblinfun_apply_left summable_on_def)
    moreover have  = (nUNIV-{0}. f C (b n *V f)) + f C (b 0 *V f)
      apply (subst infsum_Diff)
      using sum1 by auto
    moreover have  = f C (b 0 *V f) + (nUNIV-{0}. f C ((1/2 gchoose n) *R cblinfun_power S n *V f))
      unfolding b_def by simp
    moreover have  = f C (b 0 *V f) + (nUNIV-{0}. of_real (1/2 gchoose n) * (f C (cblinfun_power S n *V f)))
      by (simp add: scaleR_cblinfun.rep_eq)
    moreover have   f C (b 0 *V f) - (nUNIV-{0}. of_real (abs (1/2 gchoose n)) * (f C f)) (is _  )
    proof -
      have *: - (complex_of_real (abs (1 / 2 gchoose x)) * (f C f))
          complex_of_real (1 / 2 gchoose x) * (f C (cblinfun_power S x *V f)) for x
        apply (rule aux)
        by (auto simp: cblinfun_power_adj norm_mult fSnf selfadjoint_def
            intro!: cinner_real cinner_hermitian_real mult_left_mono Reals_mult mult_nonneg_nonneg)
      show ?thesis
        apply (subst diff_conv_add_uminus) apply (rule add_left_mono)
        apply (subst infsum_uminus[symmetric]) apply (rule infsum_mono_complex)
        apply (rule sum2)
        apply (rule sum3)
        by (rule *)
    qed
    moreover have  = f C (b 0 *V f) - (nUNIV-{0}. of_real (abs (1/2 gchoose n))) * (f C f)
      by (simp add: infsum_cmult_left')
    moreover have  = of_real (1 - (nUNIV-{0}. (abs (1/2 gchoose n)))) * (f C f)
      by (simp add: b_def left_diff_distrib infsum_of_real)
    moreover have   0 * (f C f) (is _  )
      apply (auto intro!: mult_nonneg_nonneg)
      using gbinomial_abs_has_sum_1[where a=1/2]
      by (auto simp add: infsumI)
    moreover have  = 0
      by simp
    ultimately show f C (B0 *V f)  0
      by force
  qed
  then have B  0
    by (simp add: B_def k_def scaleR_nonneg_nonneg)
  then have B = B*
    by (simp add: positive_hermitianI)
  have B0 oCL B0 = id_cblinfun + S
  proof (rule cblinfun_cinner_eqI)
    fix ψ
    define s bb where s = ψ C ((B0 oCL B0) *V ψ) and bb k = (nk. (b n *V ψ) C (b (k - n) *V ψ)) for k

    have bb k = (nk. of_real ((1 / 2 gchoose (k - n)) * (1 / 2 gchoose n)) * (ψ C (cblinfun_power S k *V ψ))) for k
      by (simp add: bb_def[abs_def] b_def cblinfun.scaleR_left cblinfun_power_adj mult.assoc
          flip: cinner_adj_right cblinfun_apply_cblinfun_compose)
    also have k = of_real (nk. ((1 / 2 gchoose n) * (1 / 2 gchoose (k - n)))) * (ψ C (cblinfun_power S k *V ψ)) for k
      apply (subst mult.commute) by (simp add: sum_distrib_right)
    also have k = of_real (1 gchoose k) * (ψ C (cblinfun_power S k *V ψ)) for k
      apply (simp only: atMost_atLeast0 gbinomial_Vandermonde)
      by simp
    also have k = of_bool (k  1) * (ψ C (cblinfun_power S k *V ψ)) for k
      by (simp add: gbinomial_1)
    finally have bb_simp: bb k = of_bool (k  1) * (ψ C (cblinfun_power S k *V ψ)) for k
      by -

    have bb_sum: bb summable_on UNIV
      apply (rule summable_on_cong_neutral[where T={..1} and g=bb, THEN iffD2])
      by (auto simp: bb_simp)

    from has_sum_b have bψ_sum: (λn. b n *V ψ) summable_on UNIV
      by (simp add: has_sum_imp_summable summable_on_cblinfun_apply_left)

    have b2_pos: (b i *V ψ) C (b j *V ψ)  0 if i0 j0 for i j
    proof -
      have gchoose_sign: (-1) ^ (i+1) * ((1/2 :: real) gchoose i)  0 if i0 for i
      proof -
        obtain j where j: Suc j = i
          using i  0 not0_implies_Suc by blast
        show ?thesis
        proof (unfold j[symmetric], induction j)
          case 0
          then show ?case
            by simp
        next
          case (Suc j)
          have (- 1) ^ (Suc (Suc j) + 1) * (1 / 2 gchoose Suc (Suc j))
               = ((- 1) ^ (Suc j + 1) * (1 / 2 gchoose Suc j)) * ((-1) * (1/2-Suc j) / (Suc (Suc j)))
            apply (simp add: gbinomial_a_Suc_n)
            by (smt (verit, ccfv_threshold) divide_divide_eq_left' divide_divide_eq_right minus_divide_right)
          also have   0
            apply (rule mult_nonneg_nonneg)
             apply (rule Suc.IH)
            apply (rule divide_nonneg_pos)
             apply (rule mult_nonpos_nonpos)
            by auto
          finally show ?case
            by -
        qed
      qed
      from S  0
      have Sn_sign: ψ C (cblinfun_power (- S) (i + j) *V ψ)  0
        by (auto intro!: cinner_pos_if_pos cblinfun_power_pos)
      have *: (- 1) ^ (i + (j + (i + j))) = (1::complex)
        by (metis Parity.ring_1_class.power_minus_even even_add power_one)

      have (b i *V ψ) C (b j *V ψ)
          = complex_of_real (1 / 2 gchoose i) * complex_of_real (1 / 2 gchoose j)
             * (ψ C (cblinfun_power S (i + j) *V ψ))
        by (simp add: b_def cblinfun.scaleR_right cblinfun.scaleR_left cblinfun_power_adj
            flip: cinner_adj_right cblinfun_apply_cblinfun_compose)
      also have  = complex_of_real ((-1)^(i+1) * (1 / 2 gchoose i)) * complex_of_real ((-1)^(j+1) * (1 / 2 gchoose j))
             * (ψ C (cblinfun_power (-S) (i + j) *V ψ))
        by (simp add: cblinfun.scaleR_left cblinfun_power_uminus * flip: power_add)
      also have   0
        apply (rule mult_nonneg_nonneg)
        apply (rule mult_nonneg_nonneg)
        using complex_of_real_nn_iff gchoose_sign that(1) apply blast
        using complex_of_real_nn_iff gchoose_sign that(2) apply blast
        by (fact Sn_sign)
      finally show ?thesis
        by -
    qed

    have s = (B0 *V ψ) C (B0 *V ψ)
      by (metis 0  B0 cblinfun_apply_cblinfun_compose cinner_adj_left positive_hermitianI s_def)
    also have  = (n. b n *V ψ) C (n. b n *V ψ)
      by (metis B0_def has_sum_b infsum_cblinfun_apply_left has_sum_imp_summable)
    also have  = (n. bb n)
      using bψ_sum bψ_sum unfolding bb_def
      apply (rule Cauchy_cinner_product_infsum[symmetric])
      using bψ_sum bψ_sum
      apply (rule Cauchy_cinner_product_summable[where X={0} and Y={0}])
      using b2_pos by auto
    also have  = bb 0 + bb 1
      apply (subst infsum_cong_neutral[where T={..1} and g=bb])
      by (auto simp: bb_simp)
    also have  = ψ C ((id_cblinfun + S) *V ψ)
      by (simp add: cblinfun_power_Suc cblinfun.add_left cinner_add_right bb_simp)
    finally show s = ψ C ((id_cblinfun + S) *V ψ)
      by -
  qed
  then have B oCL B = norm A *C (id_cblinfun + S)
    apply (simp add: k_def B_def power2_eq_square scaleR_scaleC)
    by (metis norm_imp_pos_and_ge of_real_power power2_eq_square real_sqrt_pow2)
  also have  = A
    by (metis (no_types, lifting) k_def S_def add.commute cancel_comm_monoid_add_class.diff_cancel diff_add_cancel norm_eq_zero of_real_1 of_real_mult right_inverse scaleC_diff_right scaleC_one scaleC_scaleC scaleR_scaleC)
  finally have B2A: B oCL B = A
    by -
  have BF_comm: B oCL F = F oCL B if A oCL F = F oCL A for F
  proof -
    have S oCL F = F oCL S
      by (simp add: S_def that[symmetric] cblinfun_compose_minus_right cblinfun_compose_minus_left 
          flip: cblinfun_compose_assoc)
    then have cblinfun_power S n oCL F = F oCL cblinfun_power S n for n
      apply (induction n)
       apply (simp_all add: cblinfun_power_Suc' cblinfun_compose_assoc)
      by (simp flip: cblinfun_compose_assoc)
    then have *: b n oCL F = F oCL b n for n
      by (simp add: b_def)
    have (n. b n) oCL F = F oCL (n. b n)
    proof -
      have [simp]: b summable_on UNIV
        using has_sum_b by (auto simp add: summable_on_def)
      have (n. b n) oCL F = (n. (b n) oCL F)
        apply (subst infsum_comm_additive[where f=λx. x oCL F, symmetric])
        by (auto simp: o_def isCont_cblinfun_compose_left)
      also have  = (n. F oCL (b n))
        by (simp add: * )
      also have  = F oCL (n. b n)
        apply (subst infsum_comm_additive[where f=λx. F oCL x, symmetric])
        by (auto simp: o_def isCont_cblinfun_compose_right)
      finally show ?thesis
        by -
    qed
    then have B0 oCL F = F oCL B0
      unfolding B0_def
      unfolding infsum_euclidean_eq[abs_def, symmetric]
      apply (transfer fixing: b F)
      by simp
    then show ?thesis
      by (auto simp: B_def)
  qed
  have B_closure: B  closure (cspan (range (cblinfun_power A)))
  proof (cases k = 0)
    case True
    then show ?thesis 
      unfolding B_def using closure_subset complex_vector.span_zero by auto
  next
    case False
    then have k  0
      by -
    from has_sum_b
    have limit: (sum b  B0) (finite_subsets_at_top UNIV)
      by (simp add: has_sum_def)
    have cblinfun_power (A /R k - id_cblinfun) n  cspan (range (cblinfun_power A)) for n
    proof (induction n)
      case 0
      then show ?case 
        by (auto intro!: complex_vector.span_base range_eqI[where x=0])
    next
      case (Suc n)
      define pow_n where pow_n = cblinfun_power (A /R k - id_cblinfun) n
      have pow_n_span: pow_n  cspan (range (cblinfun_power A))
        using Suc by (simp add: pow_n_def)
      have A_pow_n_span: A oCL pow_n  cspan (range (cblinfun_power A))
      proof -
        from pow_n_span
        obtain F r where finite F and F_A: F  range (cblinfun_power A)
          and pow_n_sum: pow_n = (aF. r a *C a)
          by (auto simp add: complex_vector.span_explicit)
        have A oCL a  range (cblinfun_power A) if a  F for a
        proof -
          from that obtain m where a = cblinfun_power A m
            using F_A by auto
          then have A oCL a = cblinfun_power A (Suc m)
            by (simp add: cblinfun_power_Suc')
          then show ?thesis
            by auto
        qed
        then have (aF. r a *C (A oCL a))  cspan (range (cblinfun_power A))
          by (meson basic_trans_rules(31) complex_vector.span_scale complex_vector.span_sum complex_vector.span_superset)
        moreover have A oCL pow_n = (aF. r a *C (A oCL a))
          by (simp add: pow_n_sum cblinfun_compose_sum_right flip: cblinfun.scaleC_left)
        ultimately show ?thesis
          by simp
      qed
      have cblinfun_power (A /R k - id_cblinfun) (Suc n) = (A oCL pow_n) /R k - pow_n
        by (simp add: cblinfun_power_Suc' cblinfun_compose_minus_left flip: pow_n_def)
      also from pow_n_span A_pow_n_span 
      have   cspan (range (cblinfun_power A))
        by (auto intro!: complex_vector.span_diff complex_vector.span_scale 
            simp: scaleR_scaleC)
      finally show ?case
        by -
    qed
    then have b_range: b n  cspan (range (cblinfun_power A)) for n
      by (simp add: b_def S_def scaleR_scaleC complex_vector.span_scale)
    have sum_bF: sum b F  cspan (range (cblinfun_power A)) if finite F for F
      using that apply induction
      using b_range complex_vector.span_add complex_vector.span_zero by auto
    have B0  closure (cspan (range (cblinfun_power A)))
      using limit apply (rule limit_in_closure)
      using sum_bF by (simp_all add: eventually_finite_subsets_at_top_weakI)
    also have  = closure ((λx. inverse (sqrt k) *R x) ` cspan (range (cblinfun_power A)))
      using k  0 by (simp add: scaleR_scaleC csubspace_scaleC_invariant)
    also have  = (λx. inverse (sqrt k) *R x) ` closure (cspan (range (cblinfun_power A)))
      by (simp add: closure_scaleR)
    finally show ?thesis
      apply (simp add: B_def image_def)
      using k  0 by force
  qed
  from B  0 B2A BF_comm B_closure
  show ?thesis
    by metis
qed


lemma wecken35hilfssatz:
  ― ‹Auxiliary lemma from citewecken35linearer
  P. is_Proj P  (F. F oCL (W - T) = (W - T) oCL F  F oCL P = P oCL F)
      (f. W f = 0  P f = f)
      (W = (2 *C P - id_cblinfun) oCL T)
  if WT_comm: W oCL T = T oCL W and W = W* and T = T* 
    and WW_TT: W oCL W = T oCL T
  for W T :: 'a::chilbert_space CL 'a
proof (rule exI, intro conjI allI impI)
  define P where P = Proj (kernel (W-T))
  show is_Proj P
    by (simp add: P_def)
  show thesis1: F oCL P = P oCL F if F oCL (W - T) = (W - T) oCL F for F
  proof -
    have 1: F oCL P = P oCL F oCL P if F oCL (W - T) = (W - T) oCL F for F
    proof (rule cblinfun_eqI)
      fix ψ
      have P *V ψ  space_as_set (kernel (W - T))
        by (metis P_def Proj_range cblinfun_apply_in_image)
      then have (W - T) *V P *V ψ = 0
        using kernel_memberD by blast
      then have (W - T) *V F *V P *V ψ = 0
        by (metis cblinfun.zero_right cblinfun_apply_cblinfun_compose that)
      then have F *V P *V ψ  space_as_set (kernel (W - T))
        using kernel_memberI by blast
      then have P *V (F *V P *V ψ) = F *V P *V ψ
        using P_def Proj_fixes_image by blast
      then show (F oCL P) *V ψ = (P oCL F oCL P) *V ψ
        by simp
    qed
    have 2: F* oCL (W - T) = (W - T) oCL F*
      by (metis T = T* W = W* adj_cblinfun_compose adj_minus that)
    have F oCL P = P oCL F oCL P and F* oCL P = P oCL F* oCL P
      using 1[OF that] 1[OF 2] by auto
    then show F oCL P = P oCL F
      by (metis P_def adj_Proj adj_cblinfun_compose cblinfun_assoc_left(1) double_adj)
  qed
  show thesis2: P *V f = f if W *V f = 0 for f
  proof -
    from that
    have 0 = (W *V f) C (W *V f)
      by simp
    also from W = W* have  = f C ((W oCL W) *V f)
      by (simp add: that)
    also from WW_TT have  = f C ((T oCL T) *V f)
      by simp
    also from T = T* have  = (T *V f) C (T *V f)
      by (metis cblinfun_apply_cblinfun_compose cinner_adj_left)
    finally have T *V f = 0
      by simp
    then have (W - T) *V f = 0
      by (simp add: cblinfun.diff_left that)
    then show P *V f = f
      using P_def Proj_fixes_image kernel_memberI by blast
  qed
  show thesis3: W = (2 *C P - id_cblinfun) oCL T
  proof -
    from WW_TT WT_comm have WT_binomial: (W - T) oCL (W + T) = 0
      by (simp add: cblinfun_compose_add_right cblinfun_compose_minus_left)
    have PWT: P oCL (W + T) = W + T
    proof (rule cblinfun_eqI)
      fix ψ
      from WT_binomial have (W + T) *V ψ  space_as_set (kernel (W-T))
        by (metis cblinfun_apply_cblinfun_compose kernel_memberI zero_cblinfun.rep_eq)
      then show (P oCL (W + T)) *V ψ = (W + T) *V ψ
        by (metis P_def Proj_idempotent Proj_range cblinfun_apply_cblinfun_compose cblinfun_fixes_range)
    qed
    from P_def have (W - T) oCL P = 0
      by (metis Proj_range thesis1 cblinfun_apply_cblinfun_compose cblinfun_apply_in_image
          cblinfun_eqI kernel_memberD zero_cblinfun.rep_eq)
    with PWT WT_comm thesis1 have 2 *C T oCL P = W + T
      by (metis (no_types, lifting) bounded_cbilinear.add_left bounded_cbilinear_cblinfun_compose cblinfun_compose_add_right cblinfun_compose_minus_left cblinfun_compose_minus_right eq_iff_diff_eq_0 scaleC_2)
    with  that(2) that(3) show ?thesis
      by (smt (verit, ccfv_threshold) P_def add_diff_cancel adj_Proj adj_cblinfun_compose adj_plus cblinfun_compose_id_right cblinfun_compose_minus_left cblinfun_compose_scaleC_left id_cblinfun_adjoint scaleC_2)
  qed
qed

lemma sqrt_op_pos[simp]: sqrt_op a  0
proof (cases a  0)
  case True
  from sqrt_op_existence[OF True]
  have *: b::'a CL 'a. b  0  b* oCL b = a
    by (metis positive_hermitianI)
  then show ?thesis
    using * by (smt (verit, ccfv_threshold) someI_ex sqrt_op_def)
next
  case False
  then show ?thesis
    by (simp add: sqrt_op_nonpos)
qed

lemma sqrt_op_square[simp]: 
  assumes a  0
  shows sqrt_op a oCL sqrt_op a = a
proof -
  from sqrt_op_existence[OF assms]
  have *: b::'a CL 'a. b  0  b* oCL b = a
    by (metis positive_hermitianI)
  have sqrt_op a oCL sqrt_op a = (sqrt_op a)* oCL sqrt_op a
    by (metis positive_hermitianI sqrt_op_pos)
  also have (sqrt_op a)* oCL sqrt_op a = a
    using * by (metis (mono_tags, lifting) someI_ex sqrt_op_def)
  finally show ?thesis
    by -
qed

lemma sqrt_op_unique:
  ― ‹Proof follows citewecken35linearer
  assumes b  0 and b* oCL b = a
  shows b = sqrt_op a
proof -
  have a  0
    using assms(2) positive_cblinfun_squareI by blast 
  from sqrt_op_existence[OF a  0]
  obtain sq where sq  0 and sq oCL sq = a and a_comm: a oCL F = F oCL a  sq oCL F = F oCL sq for F
    by metis
  have eq_sq: b = sq if b  0 and b* oCL b = a for b
  proof -
    have b oCL a = a oCL b
      by (metis cblinfun_assoc_left(1) positive_hermitianI that(1) that(2))
    then have b_sqrt_comm: b oCL sq = sq oCL b
      using a_comm by force
    from b  0 have b = b*
      by (simp add: assms(1) positive_hermitianI)
    have sqrt_adj: sq = sq*
      by (simp add: 0  sq positive_hermitianI)
    have bb_sqrt: b oCL b = sq oCL sq
      using b = b* sq oCL sq = a that(2) by fastforce

    from wecken35hilfssatz[OF b_sqrt_comm b = b* sqrt_adj bb_sqrt]
    obtain P where is_Proj P and b_P_sq: b = (2 *C P - id_cblinfun) oCL sq
      and Pcomm: F oCL (b - sq) = (b - sq) oCL F  F oCL P = P oCL F for F
      by metis

    have 1: sandwich (id_cblinfun - P) b = (id_cblinfun - P) oCL b
      by (smt (verit, del_insts) Pcomm is_Proj P b_sqrt_comm cblinfun_assoc_left(1) cblinfun_compose_id_left cblinfun_compose_id_right cblinfun_compose_minus_left cblinfun_compose_minus_right cblinfun_compose_zero_left diff_0_right is_Proj_algebraic is_Proj_complement is_Proj_idempotent sandwich_apply)
    also have 2:  = - (id_cblinfun - P) oCL sq
      apply (simp add: b_P_sq)
      by (smt (verit, del_insts) 0  sq is_Proj P add_diff_cancel_left' cancel_comm_monoid_add_class.diff_cancel cblinfun_compose_assoc cblinfun_compose_id_right cblinfun_compose_minus_right diff_diff_eq2 is_Proj_algebraic is_Proj_complement minus_diff_eq scaleC_2)
    also have  = - sandwich (id_cblinfun - P) sq
      by (metis (id_cblinfun - P) oCL b = - (id_cblinfun - P) oCL sq calculation cblinfun_compose_uminus_left sandwich_apply)
    also have   0
      by (simp add: 0  sq sandwich_pos)
    finally have sandwich (id_cblinfun - P) b  0
      by -
    moreover from b  0 have sandwich (id_cblinfun - P) b  0
      by (simp add: sandwich_pos)
    ultimately have sandwich (id_cblinfun - P) b = 0
      by auto
    with 1 2 have (id_cblinfun - P) oCL sq = 0
      by (metis add.inverse_neutral cblinfun_compose_uminus_left minus_diff_eq)
    with b_P_sq show b = sq
      by (metis (no_types, lifting) add.inverse_neutral add_diff_cancel_right' adj_cblinfun_compose cblinfun_compose_id_right cblinfun_compose_minus_left diff_0 diff_eq_diff_eq id_cblinfun_adjoint scaleC_2 sqrt_adj)
  qed
  
  from eq_sq have sqrt_op a = sq
    by (simp add: 0  a comparable_hermitean[unfolded selfadjoint_def])
  moreover from eq_sq have b = sq
    by (simp add: assms(1) assms(2))
  ultimately show b = sqrt_op a
    by simp
qed

lemma sqrt_op_in_closure: sqrt_op a  closure (cspan (range (cblinfun_power a)))
proof (cases a  0)
  case True
  from sqrt_op_existence[OF True]
  obtain B :: 'a CL 'a where B  0 and B oCL B = a
    and B_closure: B  closure (cspan (range (cblinfun_power a)))
    by metis
  then have sqrt_op a = B
    by (metis positive_hermitianI sqrt_op_unique)
  with B_closure show ?thesis
    by simp
next
  case False
  then have sqrt_op a = 0
    by (simp add: sqrt_op_nonpos)
  also have 0  closure (cspan (range (cblinfun_power a)))
    using closure_subset complex_vector.span_zero by blast
  finally show ?thesis
    by -
qed


lemma sqrt_op_commute:
  assumes A  0
  assumes A oCL F = F oCL A
  shows sqrt_op A oCL F = F oCL sqrt_op A
  by (metis assms(1) assms(2) positive_hermitianI sqrt_op_existence sqrt_op_unique)

lemma sqrt_op_0[simp]: sqrt_op 0 = 0
  apply (rule sqrt_op_unique[symmetric])
  by auto

lemma sqrt_op_scaleC: 
  assumes c  0 and a  0
  shows sqrt_op (c *C a) = sqrt c *C sqrt_op a
  apply (rule sqrt_op_unique[symmetric])
  using assms apply (auto simp: split_scaleC_pos_le positive_hermitianI)
  by (metis of_real_power power2_eq_square real_sqrt_pow2)

definition abs_op :: 'a::chilbert_space CL 'b::complex_inner  'a CL 'a where abs_op a = sqrt_op (a* oCL a)

lemma abs_op_pos[simp]: abs_op a  0
  by (simp add: abs_op_def positive_cblinfun_squareI sqrt_op_pos)

lemma abs_op_0[simp]: abs_op 0 = 0
  unfolding abs_op_def by auto

lemma abs_op_idem[simp]: abs_op (abs_op a) = abs_op a
  by (metis abs_op_def abs_op_pos sqrt_op_unique)

lemma abs_op_uminus[simp]: abs_op (- a) = abs_op a
  by (simp add: abs_op_def adj_uminus bounded_cbilinear.minus_left 
      bounded_cbilinear.minus_right bounded_cbilinear_cblinfun_compose)

lemma selfbutter_pos[simp]: selfbutter x  0
  by (metis butterfly_def double_adj positive_cblinfun_squareI)


lemma abs_op_butterfly[simp]: abs_op (butterfly x y) = (norm x / norm y) *R selfbutter y for x :: 'a::chilbert_space and y :: 'b::chilbert_space
proof (cases y=0)
  case False
  have abs_op (butterfly x y) = sqrt_op (cinner x x *C selfbutter y)
    unfolding abs_op_def by simp
  also have  = (norm x / norm y) *R selfbutter y
    apply (rule sqrt_op_unique[symmetric])
    using False by (auto intro!: scaleC_nonneg_nonneg simp: scaleR_scaleC power2_eq_square simp flip: power2_norm_eq_cinner)
  finally show ?thesis
    by -
next
  case True
  then show ?thesis
    by simp
qed

lemma abs_op_nondegenerate: a = 0 if abs_op a = 0
proof -
  from that
  have sqrt_op (a* oCL a) = 0
    by (simp add: abs_op_def)
  then have 0* oCL 0 = (a* oCL a)
    by (metis cblinfun_compose_zero_right positive_cblinfun_squareI sqrt_op_square)
  then show a = 0
    apply (rule_tac op_square_nondegenerate)
    by simp
qed

lemma abs_op_scaleC: abs_op (c *C a) = ¦c¦ *C abs_op a
proof -
  define aa where aa = a* oCL a
  have abs_op (c *C a) = sqrt_op (¦c¦2 *C aa)
    by (simp add: abs_op_def x_cnj_x aa_def)
  also have  = ¦c¦ *C sqrt_op aa
    by (smt (verit, best) aa_def abs_complex_def abs_nn cblinfun_compose_scaleC_left cblinfun_compose_scaleC_right complex_cnj_complex_of_real o_apply positive_cblinfun_squareI power2_eq_square scaleC_adj scaleC_nonneg_nonneg scaleC_scaleC sqrt_op_pos sqrt_op_square sqrt_op_unique)
  also have  = ¦c¦ *C abs_op a
    by (simp add: aa_def abs_op_def)
  finally show ?thesis
    by -
qed


lemma kernel_abs_op[simp]: kernel (abs_op a) = kernel a
proof (rule ccsubspace_eqI)
  fix x
  have x  space_as_set (kernel (abs_op a))  abs_op a x = 0
    using kernel_memberD kernel_memberI by blast
  also have   abs_op a x C abs_op a x = 0
    by simp
  also have   x C ((abs_op a)* oCL abs_op a) x = 0
    by (simp add: cinner_adj_right)
  also have   x C (a* oCL a) x = 0
    by (simp add: abs_op_def positive_cblinfun_squareI positive_hermitianI)
  also have   a x C a x = 0
    by (simp add: cinner_adj_right)
  also have   a x = 0
    by simp
  also have   x  space_as_set (kernel a)
    using kernel_memberD kernel_memberI by auto
  finally show x  space_as_set (kernel (abs_op a))  x  space_as_set (kernel a)
    by -
qed

definition polar_decomposition where
  ― ‹citeconway00operator, 3.9 Polar Decomposition›
  polar_decomposition A = cblinfun_extension (range (abs_op A)) (λψ. A *V inv (abs_op A) ψ) oCL Proj (abs_op A *S top)
    for A :: 'a::chilbert_space CL 'b::complex_inner

lemma 
  fixes A :: 'a :: chilbert_space CL 'b :: chilbert_space
  ― ‹citeconway00operator, 3.9 Polar Decomposition›
  shows polar_decomposition_correct: polar_decomposition A oCL abs_op A = A
    and polar_decomposition_final_space: polar_decomposition A *S top = A *S top (* Should be more precise: range (polar_decomposition A) = closure (range A) *)
    and polar_decomposition_initial_space[simp]: kernel (polar_decomposition A) = kernel A
    and polar_decomposition_partial_isometry[simp]: partial_isometry (polar_decomposition A)
proof -
  have abs_A_norm: norm (abs_op A h) = norm (A h) for h
  proof -
    have complex_of_real ((norm (A h))2) = A h C A h
      by (simp add: cdot_square_norm)
    also have  = (A* oCL A) h C h
      by (simp add: cinner_adj_left)
    also have  = ((abs_op A)* oCL abs_op A) h C h
      by (simp add: abs_op_def positive_cblinfun_squareI positive_hermitianI)
    also have  = abs_op A h C abs_op A h
      by (simp add: cinner_adj_left)
    also have  = complex_of_real ((norm (abs_op A h))2)
      using cnorm_eq_square by blast
    finally show ?thesis
      by (simp add: cdot_square_norm cnorm_eq)
  qed

  define W W' P
    where W = (λψ. A *V inv (abs_op A) ψ) 
    and W' = cblinfun_extension (range (abs_op A)) W
    and P = Proj (abs_op A *S top)

  have pdA: polar_decomposition A = W' oCL P
    by (auto simp: polar_decomposition_def W'_def W_def P_def) 

  have AA_norm: norm (W ψ) = norm ψ if ψ  range (abs_op A) for ψ
  proof -
    define h where h = inv (abs_op A) ψ
    from that have absA_h: abs_op A h = ψ
      by (simp add: f_inv_into_f h_def)
    have complex_of_real ((norm (W ψ))2) = complex_of_real ((norm (A h))2)
      using W_def h_def by blast 
    also have  = A h C A h
      by (simp add: cdot_square_norm)
    also have  = (A* oCL A) h C h
      by (simp add: cinner_adj_left)
    also have  = ((abs_op A)* oCL abs_op A) h C h
      by (simp add: abs_op_def positive_cblinfun_squareI positive_hermitianI)
    also have  = abs_op A h C abs_op A h
      by (simp add: cinner_adj_left)
    also have  = complex_of_real ((norm (abs_op A h))2)
      using cnorm_eq_square by blast
    also have  = complex_of_real ((norm ψ)2)
      using absA_h by fastforce
    finally show norm (W ψ) = norm ψ
      by (simp add: cdot_square_norm cnorm_eq)
  qed
  then have AA_norm': norm (W ψ)  1 * norm ψ if ψ  range (abs_op A) for ψ
    using that by simp

  have W_absA: W (abs_op A h) = A h for h
  proof -
    have A h = A h' if abs_op A h = abs_op A h' for h h'
    proof -
      from that have norm (abs_op A (h - h')) = 0
        by (simp add: cblinfun.diff_right)
      with AA_norm have norm (A (h - h')) = 0
        by (simp add: abs_A_norm)
      then show A h = A h'
        by (simp add: cblinfun.diff_right)
    qed
    then show ?thesis
      by (metis W_def f_inv_into_f rangeI)
  qed

  have range_subspace: csubspace (range (abs_op A))
    by (auto intro!: range_is_csubspace)

  have exP: P. is_Proj P  range ((*V) P) = closure (range ((*V) (abs_op A)))
    apply (rule exI[of _ Proj (abs_op A *S )])
    by (metis (no_types, opaque_lifting) Proj_is_Proj Proj_range Proj_range_closed cblinfun_image.rep_eq closure_closed space_as_set_top)

  have add: W (x + y) = W x + W y if x_in: x  range (abs_op A) and y_in: y  range (abs_op A) for x y
  proof -
    obtain x' y' where x = abs_op A x' and y = abs_op A y'
      using x_in y_in by blast
    then show ?thesis
      by (simp flip: cblinfun.add_right add: W_absA)
  qed

  have scale: W (c *C x) = c *C W x if x_in: x  range (abs_op A) for c x
  proof -
    obtain x' where x = abs_op A x'
      using x_in by blast
    then show ?thesis
      by (simp flip: cblinfun.scaleC_right add: W_absA)
  qed

  have cblinfun_extension_exists (range (abs_op A)) W
    using range_subspace exP add scale AA_norm'
    by (rule cblinfun_extension_exists_proj)

  then have W'_apply: W' *V ψ = W ψ if ψ  range (abs_op A) for ψ
    by (simp add: W'_def cblinfun_extension_apply that)

   have norm (W' ψ) - norm ψ = 0 if ψ  range (abs_op A) for ψ
    by (simp add: W'_apply AA_norm that)
  then have norm (W' ψ) - norm ψ = 0 if ψ  closure (range (abs_op A)) for ψ
    apply (rule_tac continuous_constant_on_closure[where S=range (abs_op A)])
    using that by (auto intro!: continuous_at_imp_continuous_on)
  then have norm_W': norm (W' ψ) = norm ψ if ψ  space_as_set (abs_op A *S top) for ψ
    using cblinfun_image.rep_eq that by force

  show correct: polar_decomposition A oCL abs_op A = A
  proof (rule cblinfun_eqI)
    fix ψ :: 'a
    have (polar_decomposition A oCL abs_op A) *V ψ = W (P (abs_op A ψ))
      by (simp add: W'_apply P_def pdA Proj_fixes_image) 
    also have  = W (abs_op A ψ)
      by (auto simp: P_def Proj_fixes_image)
    also have  = A ψ
      by (simp add: W_absA)
(*     also have ‹… = A (inv (abs_op A) (abs_op A ψ))›
      using W_def by fastforce
    also have ‹… = A ψ›
       *)
    finally show (polar_decomposition A oCL abs_op A) *V ψ = A *V ψ
      by -
  qed

  show polar_decomposition A *S top = A *S top
  proof (rule antisym)
    have *: A *S top = polar_decomposition A *S abs_op A *S top
      by (simp add: cblinfun_assoc_left(2) correct)
    also have   polar_decomposition A *S top
      by (simp add: cblinfun_image_mono)
    finally show A *S top  polar_decomposition A *S top
      by -

    have W' ψ  range A if ψ  range (abs_op A) for ψ
      using W'_apply W_def that by blast
    then have W' ψ  closure (range A) if ψ  closure (range (abs_op A)) for ψ
      using * 
      by (metis (mono_tags, lifting) P_def Proj_range Proj_fixes_image cblinfun_apply_cblinfun_compose cblinfun_apply_in_image cblinfun_compose_image cblinfun_image.rep_eq pdA that top_ccsubspace.rep_eq)
    then have W' ψ  space_as_set (A *S top) if ψ  space_as_set (abs_op A *S top) for ψ
      by (metis cblinfun_image.rep_eq that top_ccsubspace.rep_eq)
    then have polar_decomposition A ψ  space_as_set (A *S top) for ψ
      by (metis P_def Proj_range cblinfun_apply_cblinfun_compose cblinfun_apply_in_image pdA)
    then show polar_decomposition A *S top  A *S top
      using * 
      by (metis (no_types, lifting) Proj_idempotent Proj_range cblinfun_compose_image dual_order.eq_iff polar_decomposition_def)
  qed

  show partial_isometry (polar_decomposition A)
    apply (rule partial_isometryI'[where V=abs_op A *S top])
    by (auto simp add: P_def Proj_fixes_image norm_W' pdA kernel_memberD)

  have kernel (polar_decomposition A) = - (abs_op A *S top)
    apply (rule partial_isometry_initial[where V=abs_op A *S top])
    by (auto simp add: P_def Proj_fixes_image norm_W' pdA kernel_memberD)
  also have  = kernel (abs_op A)
    by (metis abs_op_pos kernel_compl_adj_range positive_hermitianI)
  also have  = kernel A
    by (simp add: kernel_abs_op)
  finally show kernel (polar_decomposition A) = kernel A
    by -
qed

lemma polar_decomposition_correct': (polar_decomposition A)* oCL A = abs_op A
  for A :: 'a :: chilbert_space CL 'b :: chilbert_space
proof -
  have polar_decomposition A* oCL A = (polar_decomposition A* oCL polar_decomposition A) oCL abs_op A
    by (simp add: cblinfun_compose_assoc polar_decomposition_correct)
  also have  = Proj (- kernel (polar_decomposition A)) oCL abs_op A
    by (simp add: partial_isometry_adj_a_o_a polar_decomposition_partial_isometry)
  also have  = Proj (- kernel A) oCL abs_op A
    by (simp add: polar_decomposition_initial_space)
  also have  = Proj (- kernel (abs_op A)) oCL abs_op A
    by simp
  also have  = Proj (abs_op A *S top) oCL abs_op A
    by (metis abs_op_pos kernel_compl_adj_range ortho_involution positive_hermitianI)
  also have  = abs_op A
    by (simp add: Proj_fixes_image cblinfun_eqI)
  finally show ?thesis
    by -
qed

lemma abs_op_adj: abs_op (a*) = sandwich (polar_decomposition a) (abs_op a)
proof -
  have pos: sandwich (polar_decomposition a) (abs_op a)  0
    by (simp add: sandwich_pos)
  have (sandwich (polar_decomposition a) (abs_op a))* oCL (sandwich (polar_decomposition a) (abs_op a))
     = polar_decomposition a oCL (abs_op a)* oCL abs_op a oCL (polar_decomposition a)*
    apply (simp add: sandwich_apply)
    by (metis (no_types, lifting) cblinfun_assoc_left(1) polar_decomposition_correct polar_decomposition_correct')
  also have  = a oCL a*
    by (metis abs_op_pos adj_cblinfun_compose cblinfun_assoc_left(1) polar_decomposition_correct positive_hermitianI)
  finally have sandwich (polar_decomposition a) (abs_op a) = sqrt_op (a oCL a*)
    using pos by (simp add: sqrt_op_unique)
  also have  = abs_op (a*)
    by (simp add: abs_op_def)
  finally show ?thesis
    by simp
qed

lemma abs_opI: 
  assumes a* oCL a = b* oCL b
  assumes a  0
  shows a = abs_op b
  by (simp add: abs_op_def assms(1) assms(2) sqrt_op_unique)

lemma abs_op_id_on_pos: a  0  abs_op a = a
  using abs_opI by force

lemma norm_abs_op[simp]: norm (abs_op a) = norm a 
  for a :: 'a::chilbert_space CL 'b::chilbert_space
proof -
  have (norm (abs_op a))2 = norm (abs_op a* oCL abs_op a)
    by simp
  also have  = norm (a* oCL a)
    by (simp add: abs_op_def positive_cblinfun_squareI positive_hermitianI)
  also have  = (norm a)2
    by simp
  finally show ?thesis
    by simp
qed


(* TODO Potentially move to Complex_Bounded_Linear_Functions and replace partial_isometry_square_proj by this.
        But the proof uses facts from this theory. *)
lemma partial_isometry_iff_square_proj:
  ― ‹@{cite conway2013course}, Exercise VIII.3.15›
  fixes A :: 'a :: chilbert_space CL 'b :: chilbert_space
  shows partial_isometry A  is_Proj (A* oCL A)
proof (rule iffI)
  show is_Proj (A* oCL A) if partial_isometry A
    by (simp add: partial_isometry_square_proj that)
next
  show partial_isometry A if is_Proj (A* oCL A)
  proof (rule partial_isometryI)
    fix h
    from that have norm (A* oCL A)  1
      using norm_is_Proj by blast
    then have normA: norm A  1 and normAadj: norm (A*)  1
      by (simp_all add: norm_AadjA abs_square_le_1)
    assume h  space_as_set (- kernel A)
    also have  = space_as_set (- kernel (A* oCL A))
      by (metis (no_types, lifting) abs_opI is_Proj_algebraic kernel_abs_op positive_cblinfun_squareI that)
    also have  = space_as_set ((A* oCL A) *S )
      by (simp add: kernel_compl_adj_range)
    finally have A* *V A *V h = h
      by (metis Proj_fixes_image Proj_on_own_range that cblinfun_apply_cblinfun_compose)
    then have norm h = norm (A* *V A *V h)
      by simp
    also have   norm (A *V h)
      by (smt (verit) normAadj mult_left_le_one_le norm_cblinfun norm_ge_zero)
    also have   norm h
      by (smt (verit) normA mult_left_le_one_le norm_cblinfun norm_ge_zero)
    ultimately show norm (A *V h) = norm h
      by simp
  qed
qed

lemma abs_op_square: (abs_op A)* oCL abs_op A = A* oCL A
  by (simp add: abs_op_def positive_cblinfun_squareI positive_hermitianI)

lemma polar_decomposition_0[simp]: polar_decomposition 0 = (0 :: 'a::chilbert_space CL 'b::chilbert_space)
proof -
  have polar_decomposition (0 :: 'a::chilbert_space CL 'b::chilbert_space) *S  = 0 *S 
    by (simp add: polar_decomposition_final_space)
  then show ?thesis
    by simp
qed

lemma polar_decomposition_unique:
  fixes A :: 'a::chilbert_space CL 'b::chilbert_space
  assumes ker: kernel X = kernel A
  assumes comp: X oCL abs_op A = A
  shows X = polar_decomposition A
proof -
  have X ψ = polar_decomposition A ψ if ψ  space_as_set (kernel A) for ψ
  proof -
    have ψ  space_as_set (kernel X)
      by (simp add: ker that)
    then have X ψ = 0
      by (simp add: kernel.rep_eq)
    moreover
    have ψ  space_as_set (kernel (polar_decomposition A))
      by (simp add: polar_decomposition_initial_space that)
    then have polar_decomposition A ψ = 0
      by (simp add: kernel.rep_eq del: polar_decomposition_initial_space)
    ultimately show ?thesis
      by simp
  qed
  then have 1: X oCL Proj (kernel A) = polar_decomposition A oCL Proj (kernel A)
    by (metis assms(1) cblinfun_compose_Proj_kernel polar_decomposition_initial_space)
  have *: abs_op A *S  = - kernel A
    by (metis (mono_tags, opaque_lifting) abs_op_pos kernel_abs_op kernel_compl_adj_range ortho_involution positive_hermitianI)
  
  have X oCL abs_op A = polar_decomposition A oCL abs_op A
    by (simp add: comp polar_decomposition_correct)
  then have X ψ = polar_decomposition A ψ if ψ  space_as_set (abs_op A *S ) for ψ
    by (simp add: cblinfun_same_on_image that)
  then have 2: X oCL Proj (- kernel A) = polar_decomposition A oCL Proj (- kernel A)
    using *
    by (metis (no_types, opaque_lifting) Proj_idempotent cblinfun_eqI lift_cblinfun_comp(4) norm_Proj_apply)
  from 1 2 have X oCL Proj (- kernel A) + X oCL Proj (kernel A)
           = polar_decomposition A oCL Proj (- kernel A) + polar_decomposition A oCL Proj (kernel A)
    by simp
  then show ?thesis
    by (simp add: Proj_ortho_compl flip: cblinfun_compose_add_right)
qed

lemma norm_cblinfun_mono:
― ‹Would logically belong in theoryComplex_Bounded_Operators.Complex_Bounded_Linear_Function
      but uses constsqrt_op from this theory in the proof.›
  fixes A B :: 'a::chilbert_space CL 'a
  assumes A  0
  assumes A  B
  shows norm A  norm B
proof -
  have B  0
    using assms by force
  have sqrtA: (sqrt_op A)* oCL sqrt_op A = A
    by (simp add: A  0 positive_hermitianI)
  have sqrtB: (sqrt_op B)* oCL sqrt_op B = B
    by (simp add: B  0 positive_hermitianI)
  have norm (sqrt_op A ψ)  norm (sqrt_op B ψ) for ψ
    apply (auto intro!: cnorm_le[THEN iffD2]
        simp: sqrtA sqrtB
        simp flip: cinner_adj_right cblinfun_apply_cblinfun_compose)
    using assms less_eq_cblinfun_def by auto
  then have norm (sqrt_op A)  norm (sqrt_op B)
    by (meson dual_order.trans norm_cblinfun norm_cblinfun_bound norm_ge_zero)
  moreover have norm A = (norm (sqrt_op A))2
    by (metis norm_AadjA sqrtA)
  moreover have norm B = (norm (sqrt_op B))2
    by (metis norm_AadjA sqrtB)
  ultimately show norm A  norm B
    by force
qed

lemma sandwich_mono: sandwich A B  sandwich A C if B  C
  by (metis cblinfun.real.diff_right diff_ge_0_iff_ge sandwich_pos that)

lemma sums_pos_cblinfun:
  fixes f :: "nat  ('b::chilbert_space CL 'b)"
  assumes f sums a
  assumes n. f n  0
  shows "a  0"
  apply (rule sums_mono_cblinfun[where f=λ_. 0 and g=f])
  using assms by auto

lemma has_sum_mono_cblinfun:
  fixes f :: "'a  ('b::chilbert_space CL 'b)"
  assumes "(f has_sum x) A" and "(g has_sum y) A"
  assumes x. x  A  f x  g x
  shows "x  y"
  using assms has_sum_mono_neutral_cblinfun by force

lemma infsum_mono_cblinfun:
  fixes f :: "'a  ('b::chilbert_space CL 'b)"
  assumes "f summable_on A" and "g summable_on A"
  assumes x. x  A  f x  g x
  shows "infsum f A  infsum g A"
  by (meson assms has_sum_infsum has_sum_mono_cblinfun)

lemma suminf_mono_cblinfun:
  fixes f :: "nat  ('b::chilbert_space CL 'b)"
  assumes "summable f" and "summable g"
  assumes x. f x  g x
  shows "suminf f  suminf g"
  using assms sums_mono_cblinfun by blast 

lemma suminf_pos_cblinfun:
  fixes f :: "nat  ('b::chilbert_space CL 'b)"
  assumes summable f
  assumes x. f x  0
  shows "suminf f  0"
  using assms sums_mono_cblinfun by blast 


lemma infsum_mono_neutral_cblinfun:
  fixes f :: "'a  ('b::chilbert_space CL 'b)"
  assumes "f summable_on A" and "g summable_on B"
  assumes x. x  AB  f x  g x
  assumes x. x  A-B  f x  0
  assumes x. x  B-A  g x  0
  shows "infsum f A  infsum g B"
  by (smt (verit, del_insts) assms(1) assms(2) assms(3) assms(4) assms(5) has_sum_infsum has_sum_mono_neutral_cblinfun)

lemma abs_op_geq: abs_op a  a if selfadjoint a
proof -
  define A P where A = abs_op a and P = Proj (kernel (A + a))
  from that have [simp]: a* = a
    by (simp add: selfadjoint_def)
  have [simp]: A  0
    by (simp add: A_def)
  then have [simp]: A* = A
    using positive_hermitianI by fastforce
  have aa_AA: a oCL a = A oCL A
    by (metis A_def A* = A abs_op_square that selfadjoint_def)
  have [simp]: P* = P
    by (simp add: P_def adj_Proj)
  have Aa_aA: A oCL a = a oCL A
    by (metis (full_types) A_def lift_cblinfun_comp(2) abs_op_def positive_cblinfun_squareI sqrt_op_commute that selfadjoint_def)

  have (A-a) ψ C (A+a) φ = 0 for φ ψ
    by (simp add: adj_minus that A* = A aa_AA Aa_aA cblinfun_compose_add_right cblinfun_compose_minus_left
        flip: cinner_adj_right cblinfun_apply_cblinfun_compose)
  then have (A-a) ψ  space_as_set (kernel (A+a)) for ψ
    by (metis A* = A adj_plus call_zero_iff cinner_adj_left kernel_memberI that selfadjoint_def)
  then have P_fix: P oCL (A-a) = (A-a)
    by (simp add: P_def Proj_fixes_image cblinfun_eqI)
  then have P oCL (A-a) oCL P = (A-a) oCL P
    by simp
  also have  = (P oCL (A-a))*
    by (simp add: adj_minus A* = A that P* = P)
  also have  = (A-a)*
    by (simp add: P_fix)
  also have  = A-a
    by (simp add: A* = A that adj_minus)
  finally have 1: P oCL (A - a) oCL P = A - a
    by -
  have 2: P oCL (A + a) oCL P = 0
    by (simp add: P_def cblinfun_compose_assoc)
  have A - a = P oCL (A - a) oCL P + P oCL (A + a) oCL P
    by (simp add: 1 2)
  also have  = sandwich P (2 *C A)
    by (simp add: sandwich_apply cblinfun_compose_minus_left cblinfun_compose_minus_right
        cblinfun_compose_add_left cblinfun_compose_add_right scaleC_2 P* = P) 
  also have   0
    by (auto intro!: sandwich_pos scaleC_nonneg_nonneg simp: less_eq_complex_def)
  finally show A  a
    by auto
qed

lemma abs_op_geq_neq: abs_op a  - a if selfadjoint a
  by (metis abs_op_geq abs_op_uminus adj_uminus that selfadjoint_def)

lemma infsum_nonneg_cblinfun:
  fixes f :: "'a  'b::chilbert_space CL 'b"
  assumes "x. x  M  0  f x"
  shows "infsum f M  0"
  apply (cases f summable_on M)
   apply (subst infsum_0_simp[symmetric])
   apply (rule infsum_mono_cblinfun)
  using assms by (auto simp: infsum_not_exists)

lemma adj_abs_op[simp]: (abs_op a)* = abs_op a
  by (simp add: positive_hermitianI) 

lemma cblinfun_image_less_eqI:
  fixes A :: 'a::complex_normed_vector CL 'b::complex_normed_vector
  assumes h. h  space_as_set S  A h  space_as_set T
  shows A *S S  T
proof -
  from assms have A ` space_as_set S  space_as_set T
    by blast
  then have closure (A ` space_as_set S)  closure (space_as_set T)
    by (rule closure_mono)
  also have  = space_as_set T
    by force
  finally show ?thesis
    apply (transfer fixing: A)
    by (simp add: cblinfun_image.rep_eq ccsubspace_leI)
qed



lemma abs_op_plus_orthogonal:
  assumes a* oCL b = 0 and a oCL b* = 0
  shows abs_op (a + b) = abs_op a + abs_op b
proof (rule abs_opI[symmetric])
  have ba: b* oCL a = 0
    apply (rule cblinfun_eqI, rule cinner_extensionality)
    apply (simp add: cinner_adj_right flip: cinner_adj_left)
    by (simp add: assms simp_a_oCL_b') 
  have abs_ab: abs_op a oCL abs_op b = 0
  proof -
    have abs_op b *S  = - kernel (abs_op b)
      by (simp add: kernel_compl_adj_range positive_hermitianI) 
    also have  = - kernel b
      by simp
    also have  = (b*) *S 
      by (simp add: kernel_compl_adj_range) 
    also have   kernel a
      apply (auto intro!: cblinfun_image_less_eqI kernel_memberI simp: )
      by (simp add: assms flip: cblinfun_apply_cblinfun_compose)
    also have  = kernel (abs_op a)
      by simp 
    finally show abs_op a oCL abs_op b = 0
      by (metis Proj_compose_cancelI cblinfun_compose_Proj_kernel cblinfun_compose_assoc cblinfun_compose_zero_left) 
  qed
  then have abs_ba: abs_op b oCL abs_op a = 0
    by (metis abs_op_pos adj_0 adj_cblinfun_compose positive_hermitianI) 
  have (a + b)* oCL (a + b) = (a*) oCL a + (b*) oCL b
    by (simp add: cblinfun_compose_add_left cblinfun_compose_add_right adj_plus assms ba)
  also have  = (abs_op a + abs_op b)* oCL (abs_op a + abs_op b)
    by (simp add: cblinfun_compose_add_left cblinfun_compose_add_right adj_plus abs_ab abs_ba flip: abs_op_square)
  finally show (abs_op a + abs_op b)* oCL (abs_op a + abs_op b) = (a + b)* oCL (a + b)
    by simp
  show 0  abs_op a + abs_op b
    by simp 
qed


definition pos_op :: 'a::chilbert_space CL 'a  'a CL 'a where
  pos_op a = (abs_op a + a) /R 2

definition neg_op :: 'a::chilbert_space CL 'a  'a CL 'a where
  neg_op a = (abs_op a - a) /R 2

lemma pos_op_pos: 
  assumes selfadjoint a
  shows pos_op a  0
  using abs_op_geq_neq[OF assms]
  apply (simp add: pos_op_def)
  by (smt (verit, best) add_le_cancel_right more_arith_simps(3) scaleR_nonneg_nonneg zero_le_divide_iff) 

lemma neg_op_pos:
  assumes selfadjoint a
  shows neg_op a  0
  using abs_op_geq[OF assms]
  by (simp add: neg_op_def scaleR_nonneg_nonneg)

lemma pos_op_neg_op_ortho:
  assumes selfadjoint a
  shows pos_op a oCL neg_op a = 0
  apply (auto intro!: simp: pos_op_def neg_op_def cblinfun_compose_add_left cblinfun_compose_minus_right)
  by (metis (no_types, opaque_lifting) Groups.add_ac(2) abs_op_def abs_op_pos abs_op_square assms cblinfun_assoc_left(1) positive_cblinfun_squareI positive_hermitianI selfadjoint_def sqrt_op_commute) 


lemma pos_op_plus_neg_op: pos_op a + neg_op a = abs_op a
  by (simp add: pos_op_def neg_op_def scaleR_diff_right scaleR_add_right pth_8)

lemma pos_op_minus_neg_op: pos_op a - neg_op a = a
  by (simp add: pos_op_def neg_op_def scaleR_diff_right scaleR_add_right pth_8)

lemma pos_op_neg_op_unique:
  assumes bca: b - c = a
  assumes b  0 and c  0
  assumes bc: b oCL c = 0
  shows b = pos_op a and c = neg_op a
proof -
  from bc have cb: c oCL b = 0
    by (metis adj_0 adj_cblinfun_compose assms(2) assms(3) positive_hermitianI) 
  from b  0 have [simp]: b* = b
    by (simp add: positive_hermitianI) 
  from c  0 have [simp]: c* = c
    by (simp add: positive_hermitianI) 
  have bc_abs: b + c = abs_op a
  proof -
    have (b + c)* oCL (b + c) = b oCL b + c oCL c
      by (simp add: cblinfun_compose_add_left cblinfun_compose_add_right bc cb adj_plus)
    also have  = (b - c)* oCL (b - c)
      by (simp add: cblinfun_compose_minus_left cblinfun_compose_minus_right bc cb adj_minus)
    also from bca have  = a* oCL a
      by blast
    finally show ?thesis
      apply (rule abs_opI)
      by (simp add: b  0 c  0) 
  qed
  from arg_cong2[OF bca bc_abs, of plus]
    arg_cong2[OF pos_op_minus_neg_op[of a] pos_op_plus_neg_op[of a], of plus, symmetric]
  show b = pos_op a
    by (simp flip: scaleR_2)
  from arg_cong2[OF bc_abs bca, of minus]
    arg_cong2[OF pos_op_plus_neg_op[of a] pos_op_minus_neg_op[of a], of minus, symmetric]
  show c = neg_op a
    by (simp flip: scaleR_2)
qed


lemma pos_imp_selfadjoint: a  0  selfadjoint a
  using positive_hermitianI selfadjoint_def by blast

lemma abs_op_one_dim: abs_op x = one_dim_iso (abs (one_dim_iso x :: complex))
  by (metis (mono_tags, lifting) abs_opI abs_op_scaleC of_complex_def one_cblinfun_adj one_comp_one_cblinfun one_dim_iso_is_of_complex one_dim_iso_of_one one_dim_iso_of_zero one_dim_loewner_order one_dim_scaleC_1 zero_less_one_class.zero_le_one)


lemma pos_selfadjoint: selfadjoint a if a  0
  using adj_0 comparable_hermitean selfadjoint_def that by blast

lemma one_dim_loewner_order_strict: A > B  one_dim_iso A > (one_dim_iso B :: complex) for A B :: 'a CL 'a::{chilbert_space, one_dim}
  by (auto simp: less_cblinfun_def one_dim_loewner_order)

lemma one_dim_cblinfun_zero_le_one: 0 < (1 :: 'a::one_dim CL 'a)
  by (simp add: one_dim_loewner_order_strict)
lemma one_dim_cblinfun_one_pos: 0  (1 :: 'a::one_dim CL 'a)
  by (simp add: one_dim_loewner_order)

lemma Proj_pos[iff]: Proj S  0
  apply (rule positive_cblinfun_squareI[where B=Proj S])
  by (simp add: adj_Proj)

lemma abs_op_Proj[simp]: abs_op (Proj S) = Proj S
  by (simp add: abs_op_id_on_pos)



lemma diagonal_operator_pos:
  assumes x. f x  0
  shows diagonal_operator f  0
proof (cases bdd_above (range (λx. cmod (f x))))
  case True
  have [simp]: csqrt (f x) = sqrt (cmod (f x)) for x
    by (simp add: Extra_Ordered_Fields.complex_of_real_cmod assms abs_pos of_real_sqrt)
  have bdd: bdd_above (range (λx. sqrt (cmod (f x))))
  proof -
    from True obtain B where cmod (f x)  B for x
      by (auto simp: bdd_above_def)
    then show ?thesis
      by (auto intro!: bdd_aboveI[where M=sqrt B] simp: )
  qed
  show ?thesis
    apply (rule positive_cblinfun_squareI[where B=diagonal_operator (λx. csqrt (f x))])
    by (simp add: assms diagonal_operator_adj diagonal_operator_comp bdd complex_of_real_cmod abs_pos
        flip: of_real_mult)
next
  case False
  then show ?thesis
    by (simp add: diagonal_operator_invalid)
qed

lemma abs_op_diagonal_operator: 
  abs_op (diagonal_operator f) = diagonal_operator (λx. abs (f x))
proof (cases bdd_above (range (λx. cmod (f x))))
  case True
  show ?thesis
    apply (rule abs_opI[symmetric])
    by (auto intro!: diagonal_operator_pos abs_nn simp: True diagonal_operator_adj diagonal_operator_comp cnj_x_x)
next
  case False
  then show ?thesis
    by (simp add: diagonal_operator_invalid)
qed



end