Theory Jacobi_Triple_Product

section ‹The Jacobi Triple Product›
theory Jacobi_Triple_Product
  imports Theta_Functions "Lambert_Series.Lambert_Series_Library"
begin

(*<*)
(* TODO: notation for the "old" deprecated infinite sums; why is this even still active *)
no_notation Infinite_Set_Sum.abs_summable_on (infix "abs'_summable'_on" 50)
(*>*)

subsection ‹Versions for Jacobi's theta function›

unbundle qpochhammer_inf_notation

text ‹
  The following follows the short proof given by Andrews~\cite{andrews},
  which makes use of two series expansions for $(a; b)_\infty$ and $1/(a; b)\infty$ 
  due to Euler.

  We prove it for Jacobi's theta function and derive a version for Ramanujan's later on.
  One could possibly also adapt the prove to work for Ramanujan's version directly, which
  makes the transfer to Jacobi's function a bit easier. However, I chose to do it for
  Jacobi's version first in order to stay closer to the text by Andrews.

  The proof is fairly straightforward; the only messy part is proving the absolute convergence
  of the double sum (which Andrews does not bother doing). This is necessary in order to
  allow the exchange of the order of summation.
›
theorem jacobi_theta_nome_triple_product_complex:
  fixes w q :: complex
  assumes "w  0" "norm q < 1"
  shows   "jacobi_theta_nome w q = (q2 ; q2) * (-q*w ; q2) * (-q/w ; q2)"
proof -
  text ‹
    We first establish the identity for real $w$ and $q$ with the somewhat arbitrarily
    chosen bounds $q\in (0,\frac{1}{4})$ and $w\in (\frac{1}{2}, 1)$. This is considerably
    more restrictive than Andrews's version, but I was not able to prove absolute convergence
    of the sum for his bounds.

    It does not matter anyway, since we will extend it to the full domain of θ› by
    analytic continuation later anyway.
  ›
  have eq_real: "jacobi_theta_nome w q = (q2; q2) * (-q*w; q2) * (-q/w; q2)"
    if wq: "0 < q" "q < 1/4" "1/2 < w" "w < 1" for q w :: real
  proof -
    have q2: "q^2 < 1"
      using wq by (simp add:  power_less_one_iff)
    define tri where "tri = (λn::nat. n*(n+1) div 2)"
    have [simp]: "2 * (n * (n - Suc 0) div 2) = n * (n - Suc 0)" for n
      by (subst dvd_mult_div_cancel) auto
    from wq have [simp]: "w  0" "q  0"
      by auto
    have [simp]: "(q2; q2)  0"
    proof
      assume "(q2; q2) = 0"
      then obtain n where "q^2 * q ^ (2*n) = 1"
        using q2 wq by (auto simp: qpochhammer_inf_zero_iff power_mult)
      hence "q ^ (2*n + 2) = 1"
        by (simp add: eval_nat_numeral)
      moreover have "q ^ (2*n + 2) < 1"
        by (subst power_less_one_iff) (use wq in auto)
      ultimately show False
        by linarith
    qed

    have "((λn. (q2)^(n*(n-1) div 2) * (-w*q)^n / (k=1..n. (q2)^k - 1)) has_sum (-w*q;q2)) UNIV"
      by (intro norm_summable_imp_has_sum sums_qpochhammer_inf_real norm_summable_qpochhammer_inf)
         (use wq q2 mult_strict_mono[of w 1 q 1] in auto)
    hence "((λn. q ^ (n*(n-1)) * (q^n * (-w)^n) / (k=1..n. q^(2*k) - 1)) has_sum (-w*q; q2)) UNIV"
      by (simp add: norm_power power_less_one_iff power_minus' power_mult_distrib mult_ac
               flip: power_mult)
    also have "(λn. q ^ (n*(n-1)) * (q^n * (-w)^n) / (k=1..n. q^(2*k) - 1)) = 
               (λn. q ^ (n^2) * w^n * (q^(2*n+2); q^2) / (q^2;q^2))" (is "?lhs = ?rhs")
    proof
      fix n :: nat
      have "(q2 ; q2) = (k<n. 1 - q^(2*k+2)) * (q^(2*n+2) ; q2)"
        using qpochhammer_inf_mult_power_q[of "q^2" "q^2" n] q2
        by (simp add: power_add power2_eq_square power_mult_distrib mult_ac mult_2_right
                      qpochhammer_nonneg_def)
      also have "(k<n. 1 - q^(2*k+2)) = (k=1..n. 1 - q^(2*k))"
        by (rule prod.reindex_bij_witness[of _ "λk. k-1" Suc]) auto
      also have "(k=1..n. 1 - q^(2*k)) = (-1)^n * (k=1..n. q^(2*k)-1)"
        by (subst prod_diff_swap) auto
      finally have "?lhs n = q ^ (n*(n-1)) * q^n * w^n * (q^(2*n+2); q^2) / (q^2;q^2)"
        by (auto simp: divide_simps power_minus')
      also have "q^(n*(n-1)) * q^n = q^(n*(n-1) + n)"
        by (simp add: power_add)
      also have "n*(n-1) + n = n^2"
        by (simp add: power2_eq_square algebra_simps)
      finally show "?lhs n = ?rhs n" .
    qed
    finally have "((λn. (q2; q2) * (q^(n^2) * w^n * (q^(2*n+2); q2) / (q2;q2)))
                     has_sum (q2; q2) * (-w*q;q2)) UNIV"
      by (rule has_sum_cmult_right)
    hence "((λn. q^(n^2) * w^n * (q^(2*n+2); q2)) has_sum (q2;q2) * (-w*q;q2)) UNIV"
      by simp
    also have "?this  ((λn. q powi (n^2) * w powi n * (q powi (2*n+2); q2)) has_sum 
                           (q2;q2) * (-w*q;q2)) {0..}"
      by (intro has_sum_reindex_bij_witness[of _ nat int])
         (auto simp: power_int_def nat_add_distrib nat_mult_distrib nat_power_eq)
    also have "  ((λn. q powi (n^2) * w powi n * (q powi (2*n+2); q2)) has_sum 
                           (q2;q2) * (-w*q;q2)) UNIV"
    proof (rule has_sum_cong_neutral)
      fix n :: int
      assume "n  UNIV - {0..}"
      hence n: "n < 0"
        by simp
      define k where "k = nat (-n-1)"
      have "q powi (2 * int k) = (q ^ 2) ^ k"
        by (auto simp: power_int_def power_mult nat_mult_distrib)
      hence "q powi (2 * n + 2) * q2 ^ k = q powi (2*(n + 1 + int k))"
        by (simp add: power_int_add flip: power_mult power_int_mult)
      also have "n + 1 + int k = 0"
        using n by (auto simp: k_def)
      finally have "k. q powi (2 * n + 2) * q2 ^ k = 1"
        by (intro exI[of _ k]) auto
      thus "q powi (n^2) * w powi n * (q powi (2 * n + 2) ; q2) = 0"
        using q2 by (auto simp: qpochhammer_inf_zero_iff intro!: exI[of _ "-1"])
    qed auto
    finally have "((λn. q powi (n^2) * w powi n * (q powi (2*n+2); q2)) has_sum 
                    (q2;q2) * (-w*q;q2)) UNIV" .

    text ‹
      Brace yourselves: now we have to prove absolute convergence of the double sum.
      We use a crude bound for the inner sum, at which point the outer sum is just
      a geometric one that obviously converges.
    ›
    have "((λ(n,m). q powi n2 * w powi n * ((-1)^m * q powi (m^2+m+2*n*m) / (i=1..m. 1 - q^(2*i)))) has_sum
            (q2; q2) * (-w*q ; q2)) (UNIV × UNIV)"
    proof (rule has_sum_SigmaI; (unfold prod.case)?)
      show "((λn. q powi (n^2) * w powi n * (q powi (2*n+2); q2)) 
              has_sum (q2;q2) * (-w*q;q2)) UNIV"
        by fact
    next
      fix n :: int
      have "((λm. q2 ^ (m*(m-1) div 2) * (q powi (2*n+2))^m / (k=1..m. q2^k-1)) has_sum
                      (q powi (2*n+2);q2)) UNIV"
        by (intro norm_summable_imp_has_sum sums_qpochhammer_inf_real
                  norm_summable_qpochhammer_inf) (use q2 in auto)
      also have "(λm. (q2) ^ (m*(m-1) div 2) * (q powi (2*n+2))^m / (k=1..m. q2^k-1)) =
                 (λm. (-1)^m * q powi (m^2+m+2*n*m) / (k=1..m. 1-q^(2*k)))" (is "?lhs = ?rhs")
      proof
        fix m :: nat
        have "(q2) ^ (m*(m-1) div 2) * (q powi (2*n+2))^m = 
              q powi (int m * int (m-1) + int m * (2*n+2))" 
          by (simp add: power_int_add power_int_nonneg_exp nat_mult_distrib power_int_power' 
                        power_mult_distrib ring_distribs mult_ac
                   del: power_Suc flip: power_mult)
        also have "int m * int (m-1) + int m * (2*n+2) = int m ^ 2 + int m + 2 * n * int m"
          by (cases m) (simp_all add: algebra_simps power2_eq_square)
        finally have "q2 ^ (m * (m - 1) div 2) * (q powi (2 * n + 2)) ^ m =
                      q powi (int m ^ 2 + int m + 2 * n * int m)" .
        moreover have "(k=1..m. q2^k-1) = (-1)^m * (k=1..m. 1 - q^(2*k))"
        proof -
          have "(k=1..m. q2^k-1) = (k=1..m. q^(2*k)-1)"
            by (simp add: power_mult)
          also have " = (-1)^m * (k=1..m. 1 - q^(2*k))"
            by (subst prod_diff_swap) auto
          finally show ?thesis .
        qed
        ultimately show "?lhs m = ?rhs m"
          by (simp add: divide_simps)
      qed
      finally show "((λm. q powi (n^2) * w powi n * ((-1)^m * q powi (m^2+m+2*n*m) / (k=1..m. 1 - q^(2*k)))) has_sum
                       (q powi (n^2) * w powi n * (q powi (2*n+2);q2))) UNIV"
        by (rule has_sum_cmult_right)
    next
      have "(λ(n, m). q powi n2 * w powi n * (q powi ((int m)2 + int m + 2 * n * int m) /
              (i = 1..m. 1 - q ^ (2 * i)))) summable_on UNIV × UNIV"
      proof (rule summable_on_comparison_test; (safe)?)
        fix n :: int and m :: nat
        have "q powi (n2) * w powi n * (q powi (m2 + m + 2*n*m) / (i=1..m. 1-q^(2*i))) 
              q powi (n2) * w powi n * (q powi (m2 + m + 2*n*m) / (i=1..m. 1-q^1))"
          by (intro mult_left_mono divide_left_mono prod_mono zero_le_power_int conjI
                       diff_left_mono power_decreasing mult_nonneg_nonneg prod_pos mult_pos_pos)
             (use wq in auto simp: power_less_one_iff power_le_one_iff)
        also have " = q powi (n2) * w powi n * (q powi (m2 + m + 2*n*m) / (1-q)^m)"
          by simp
        also have " = q powi (n+m)^2 * w powi n * (q / (1-q))^m"
          by (simp add: power2_eq_square algebra_simps power_int_add power_divide)
        also have " = q powi (n+m)^2 * w powi (n+m) * (q / (w * (1 - q)))^m"
          using wq by (simp add: power_int_add divide_simps mult_ac)
        also have "  w powi (n+m)^2 * w powi (n+m) * (q / (w * (1 - q)))^m"
          by (rule mult_right_mono power_int_mono)+ (use wq in auto)
        also have " = w powi ((n+m)^2+(n+m)) * (q/(w*(1-q)))^m"
          by (simp add: power2_eq_square power_int_add)             
        finally show "q powi n2 * w powi n * (q powi (int m ^ 2 + int m + 2*n * int m) / 
                         (i=1..m. 1-q^(2*i)))
                       (case (n, m) of (n, m)  w powi ((n+m)^2+(n+m)) * (q/(w*(1-q)))^m)"
          by simp
      next
        have "(λ(n,m). w powi (n^2+n) * (q/(w*(1-q)))^m) summable_on UNIV × UNIV"
        proof (rule summable_on_SigmaI; (unfold prod.case)?)
          fix n :: int
          have "q < (1 / 2) * (3 / 4)"
            using wq by simp
          also have "  w * (1 - q)"
            by (intro mult_mono) (use wq in auto)
          finally have "(λm. (q/(w*(1-q)))^m) sums (1 / (1 - q/(w*(1-q))))"
            by (intro geometric_sums) (use wq in auto simp: power_le_one_iff power_less_one_iff)
          hence "((λm. (q/(w*(1-q)))^m) has_sum (1 / (1 - q/(w*(1-q))))) UNIV"
            by (intro sums_nonneg_imp_has_sum) (use wq in auto)
          thus "((λm. w powi (n2 + n) * (q / (w * (1 - q))) ^ m) 
                  has_sum (w powi (n2 + n)) * (1/(1-q/(w*(1-q))))) UNIV"
            by (intro has_sum_cmult_right)
        next
          have "summable (λn. w ^ n)"
            by (intro summable_mult2 summable_geometric) (use wq in auto)
          hence "(λn. w ^ n) summable_on UNIV"
            by (rule summable_nonneg_imp_summable_on) (use wq in auto)
          hence "(λn. w ^ n * (1 / (1 - q / (w * (1 - q))))) summable_on UNIV"
            by (intro summable_on_cmult_left)
          hence "(λn. w ^ n * (1 / (1 - q / (w * (1 - q))))) summable_on range (λn. n*(n+1))"
            by (rule summable_on_subset) auto
          also have "?this  (λn. w ^ (n*(n+1)) * (1 / (1 - q/(w*(1-q))))) summable_on UNIV"
          proof (subst summable_on_reindex)
            have "strict_mono (λn::nat. n * (n + 1))"
              by (intro strict_monoI_Suc) auto
            thus "inj (λn::nat. n * (n + 1))"
              by (rule strict_mono_imp_inj_on)
          qed (simp_all only: o_def)
          also have "  (λn. w powi (n*(n+1)) * (1 / (1 - q/(w*(1-q))))) summable_on {0..}"
            by (rule summable_on_reindex_bij_witness[of _ nat int])
               (auto simp: power_int_def nat_mult_distrib nat_add_distrib)
          finally have summable1: "(λn. w powi (n*(n+1)) * (1/(1-q/(w*(1-q))))) summable_on {0..}" .
          also have "?this  (λn. w powi (n*(n+1)) * (1/(1-q/(w*(1-q))))) summable_on {..-1}"
            by (rule summable_on_reindex_bij_witness[of _ "λn. -n-1" "λn. -n-1"]) 
               (auto simp: algebra_simps)
          finally have summable2: "(λn. w powi (n*(n+1)) * (1/(1-q/(w*(1-q))))) summable_on {..-1}" .
          have "(λn. w powi (n*(n+1)) * (1/(1-q/(w*(1-q))))) summable_on ({..-1}  {0..})"
            by (intro summable_on_union summable1 summable2)
          also have "{..-1}  {0::int..} = UNIV"
            by auto
          finally show "(λn. w powi (n^2+n) * (1/(1-q/(w*(1-q))))) summable_on UNIV"
            by (simp add: power2_eq_square algebra_simps)
        qed (use wq in auto)
        also have "?this  ((λ(n,m). w powi ((n + int m)2 + (n + int m)) * (q/(w*(1-q))) ^ m)
                                summable_on UNIV × UNIV)"
          by (rule summable_on_reindex_bij_witness
                [of _ "λ(n,m). (n + int m, m)" "λ(n,m). (n - int m, m)"]) auto
        finally show "(λ(n,m). w powi ((n+m)2 + (n+m)) * (q/(w*(1-q)))^m) summable_on UNIV × UNIV" .
      qed (use wq in auto intro!: mult_nonneg_nonneg divide_nonneg_pos prod_pos
                           simp: power_less_one_iff)
      hence "(λ(n, m). q powi n2 * w powi n * ((-1) ^ m * q powi ((int m)2 + int m + 2 * n * int m) /
              (i=1..m. 1 - q ^ (2 * i)))) abs_summable_on UNIV × UNIV"
        using wq by (simp add: case_prod_unfold abs_mult power_abs
                               power_int_abs abs_prod power_le_one_iff)
      thus "(λ(n, m). q powi n2 * w powi n * ((- 1) ^ m * q powi ((int m)2 + int m + 2 * n * int m) /
              (i = 1..m. 1 - q ^ (2 * i)))) summable_on UNIV × UNIV"
        by (rule abs_summable_summable)
    qed

    also have "(λ(n,m). q powi n2 * w powi n * ((-1)^m * q powi (int m ^ 2 + int m + 2 * n * int m) /
                   (i=1..m. 1 - q^(2*i)))) =
               (λ(n,m). ((-1)^m * (q/w)^m / (i=1..m. 1 - q^(2*i))) *
                        (q powi ((n+m)^2) * w powi (n+m)))"
      by (simp add: power2_eq_square field_simps power_int_add)
    also have "( has_sum ((q2; q2) * (-w*q; q2))) (UNIV × UNIV) 
                 ((λ(n,m). ((-1)^m * (q/w)^m / (i=1..m. 1 - q^(2*i))) * (q powi (n^2) * w powi n))
                    has_sum ((q2; q2) * (-w*q; q2))) (UNIV × UNIV)"
      by (rule has_sum_reindex_bij_witness[of _ "λ(n, m). (n - int m, m)" "λ(n, m). (n + int m, m)"])
         auto
    finally have sum:
      "((λ(n,m). ((-q/w)^m / (i=1..m. 1 - q^(2*i))) * (q powi (n^2) * w powi n))
         has_sum ((q2; q2) * (-w*q; q2))) (UNIV × UNIV)" 
      by (simp add: power_minus')

    have "((λn. inverse (-q/w; q2) * (q powi (n^2) * w powi n)) 
            has_sum (q2 ; q2) * (- w * q ; q2)) UNIV"
    proof (rule has_sum_SigmaD[OF sum]; unfold prod.case)
      fix n :: int
      have "((λm. (-q/w)^m / (i=1..m. 1 - (q2)^i)) has_sum (inverse (-q/w; q2))) UNIV"
        by (intro norm_summable_imp_has_sum sums_inverse_qpochhammer_inf_real
                  norm_summable_inverse_qpochhammer_inf)
           (use q2 wq in auto simp: norm_divide)
      also have "(λi. (q^2)^i) = (λi. q^(2*i))"
        by (simp add: power_mult)
      finally show "((λm. ((-q/w)^m / (i=1..m. 1 - q^(2*i))) * (q powi (n^2) * w powi n)) 
                       has_sum (inverse (-q/w; q2) * (q powi (n^2) * w powi n))) UNIV"
        by (rule has_sum_cmult_left)
    qed
    hence "((λn. (-q/w; q2) * (inverse (-q/w; q2) * (q powi (n^2) * w powi n)))
              has_sum (-q/w;q2) * ((q2;q2) * (-w*q;q2))) UNIV"
      by (rule has_sum_cmult_right)
    moreover have "(-q/w; q2)  0"
    proof
      assume "(-q/w; q2) = 0"
      then obtain n where "q ^ (2*n+1) / w = -1"
        using q2 by (auto simp: qpochhammer_inf_zero_iff power_mult)
      have "norm (q ^ (2 * n + 1) / w) = q ^ (2*n+1) / w"
        using wq by (simp add: norm_divide norm_power)
      also have "  q ^ 1 / w"
        by (intro divide_right_mono power_decreasing) (use wq in auto)
      also have " < 1"
        using wq by simp
      also note q ^ (2*n+1) / w = -1
      finally show False
        by simp
    qed
    ultimately have "((λn. w powi n * q powi (n^2)) has_sum 
                       (q2; q2) * (-q*w; q2) * (-q/w; q2)) UNIV"
      by (simp add: divide_simps mult_ac)
    moreover have "((λn. w powi n * q powi (n^2)) has_sum jacobi_theta_nome w q) UNIV"
      by (rule has_sum_jacobi_theta_nome) (use wq in auto)
    ultimately show "jacobi_theta_nome w q = (q2; q2) * (-q*w; q2) * (-q/w; q2)"
      using has_sum_unique by blast
  qed

  text ‹
    We perform analytic continuation to extend the identity to all $w$:
  ›
  have eq_real2:
     "jacobi_theta_nome w (complex_of_real q) = 
        (of_real q ^ 2; of_real q ^ 2) * (-of_real q * w; of_real q ^ 2) * 
        (-of_real q / w; of_real q ^ 2)" 
    if wq: "w  0" "q  {0<..<1/4}" for w :: complex and q :: real
  proof -
    define r where "r = (2/3 :: real)"
    have r: "q < r" "1/2 < r" "r < 1"
      using wq by (simp_all add: r_def)
    define f where
      "f = (λw. jacobi_theta_nome w (complex_of_real q) - 
                (of_real q ^ 2; of_real q ^ 2) * (-of_real q * w; of_real q ^ 2) * 
                (-of_real q / w; of_real q ^ 2))"
    have "f w = 0"
    proof (rule analytic_continuation[of f])
      show "f holomorphic_on (-{0})"
        unfolding f_def using wq
        by (auto intro!: holomorphic_intros simp: norm_power power_less_one_iff)
    next
      show "of_real ` {1/2<..<1}  -{0}"
        using wq by auto
    next
      show "of_real r islimpt complex_of_real ` {1/2<..<1}" using r
        by (intro islimpt_isCont_image continuous_intros)
           (auto simp: eventually_at_filter open_imp_islimpt simp del: of_real_add)
    next
      show "f w = 0" if "w  complex_of_real ` {1/2<..<1}" for w :: complex
      proof -
        from that obtain x where x: "w = complex_of_real x" "x  {1/2<..<1}"
          by auto
        have "0 = complex_of_real (jacobi_theta_nome x q - (q2 ; q2) * (-q*x; q2) * (-q/x ; q2))"
          by (subst eq_real) (use wq x in auto)
        also have " = f w" using x wq
          by (simp add: f_def power_less_one_iff 
                   flip: jacobi_theta_nome_of_real qpochhammer_inf_of_real)
        finally show ?thesis ..
      qed
    qed (use wq r in auto simp: connected_punctured_universe)
    thus ?thesis
      by (simp add: f_def)
  qed

  text ‹
    And another analytic continuation to extend it to all $q$:
  ›
  show "jacobi_theta_nome w q = (q ^ 2; q ^ 2) * (-q * w; q ^ 2) * (-q / w; q ^ 2)" 
  proof -
    note wq = assms
    define f where
      "f = (λq. jacobi_theta_nome w q - 
                (q ^ 2; q ^ 2) * (-q * w; q ^ 2) * (-q / w; q ^ 2))"
    have "f q = 0"
    proof (rule analytic_continuation[of f])
      show "f holomorphic_on (ball 0 1)"
        unfolding f_def using wq
        by (auto intro!: holomorphic_intros simp: norm_power power_less_one_iff)
    next
      show "of_real ` {0<..<1/4}  ball (0::complex) 1"
        using wq by auto
    next
      show "of_real (1/8) islimpt complex_of_real ` {0<..<1/4}"
        by (intro islimpt_isCont_image continuous_intros)
           (auto simp: eventually_at_filter open_imp_islimpt complex_eq_iff)
    next
      show "f q = 0" if "q  complex_of_real ` {0<..<1/4}" for q :: complex
      proof -
        from that obtain x where x: "q = complex_of_real x" "x  {0<..<1/4}"
          by auto
        have "0 = jacobi_theta_nome w x - (of_real x ^ 2 ; of_real x ^ 2) *
                  (-of_real x * w; of_real x ^ 2) * (-of_real x / w ; of_real x ^ 2)"
          by (subst eq_real2) (use wq x in auto)
        also have " = f q" using x wq
          by (simp add: f_def power_less_one_iff 
                   flip: jacobi_theta_nome_of_real qpochhammer_inf_of_real)
        finally show ?thesis ..
      qed
    qed (use wq in auto)
    thus ?thesis
      by (simp add: f_def)
  qed
qed

lemma jacobi_theta_nome_triple_product_real:
  fixes w q :: real
  assumes "w  0" "¦q¦ < 1"
  shows   "jacobi_theta_nome w q = (q2 ; q2) * (-q*w ; q2) * (-q/w ; q2)"
proof -
  define q' w' where "q' = complex_of_real q" and "w' = complex_of_real w"
  have "complex_of_real (jacobi_theta_nome w q) = jacobi_theta_nome w' q'"
    by (simp add: jacobi_theta_nome_of_real w'_def q'_def)
  also have " = (q'2 ; q'2) * (-q'*w' ; q'2) * (-q'/w' ; q'2)"
    by (rule jacobi_theta_nome_triple_product_complex)
       (use assms in simp_all add: q'_def w'_def)
  also have " = of_real ((q2 ; q2) * (-q*w ; q2) * (-q/w ; q2))" using assms 
    by (simp add: q'_def w'_def power_less_one_iff abs_square_less_1 flip: qpochhammer_inf_of_real)
  finally show ?thesis
    by (simp only: of_real_eq_iff)
qed


subsection ‹Version of Ramanujan's theta function›

text ‹
  The triple product for Ramanujan's theta function, which follows easily from the above one,
  has a particularly elegant form:
  \[f(a,b) = (-a; ab)_\infty\,(-b; ab)_\infty\, (ab; ab)_\infty\]
  It follows directly from the version for Jacobi's theta function, although I again have to
  employ analytic continuation to avoid dealing with the branch cuts when converting 
  Ramanujan's theta function to Jacobi's.
›
theorem ramanujan_theta_triple_product_complex:
  fixes a b :: complex
  assumes "norm (a * b) < 1"
  shows   "ramanujan_theta a b = (-a ; a*b) * (-b ; a*b) * (a*b ; a*b)"
proof -
  have real_eq1: "ramanujan_theta a b = (-a ; a * b) * (-b ; a * b) * (a * b ; a * b)"
    if ab: "a > 0" "b > 0" "a * b < 1" for a b :: real
  proof -
    define w q where "w = sqrt (a/b)" and "q = sqrt (a*b)"
    have [simp]: "w  0"
      using ab by (auto simp: w_def sgn_if)
    have q: "¦q¦ < 1"
      using ab by (simp_all add: q_def abs_mult flip: real_sqrt_abs')
  
    have "ramanujan_theta a b = jacobi_theta_nome w q"
      using ab by (auto simp: sgn_if jacobi_theta_nome_def w_def q_def real_sqrt_mult real_sqrt_divide)
    also have " = (q2 ; q2) * (- q * w ; q2) * (- q / w ; q2)"
      by (rule jacobi_theta_nome_triple_product_real) (use q in auto)
    also have " = (-a ; a * b) * (-b ; a * b) * (a * b ; a * b)"
      using ab by (simp add: mult_ac q_def w_def real_sqrt_mult real_sqrt_divide power_mult_distrib)
    finally show ?thesis .
  qed

  have real_eq2: "ramanujan_theta a b = (-a ; a * b) * (-b ; a * b) * (a * b ; a * b)"
    if ab: "a > 0" "norm (a * b) < 1" for a :: real and b :: complex
  proof -
    define f :: "complex  complex"
      where "f = (λb. ramanujan_theta (of_real a) b - (-of_real a ; of_real a * b) * 
                      (-b ; of_real a * b) * (of_real a * b ; of_real a * b))"

    have "f b = 0"
    proof (rule analytic_continuation[where f = f])
      show "f holomorphic_on (ball 0 (1 / a))"
        unfolding f_def using ab by (intro holomorphic_intros) (auto simp: norm_mult field_simps)
    next
      show "complex_of_real ` {0<..<1/a}  ball 0 (1/a)"
       and "complex_of_real (1/(2*a))  ball 0 (1/a)"
        using ab by (auto simp: norm_mult norm_divide field_simps)
    next
      show "complex_of_real (1 / (2 * a)) islimpt
            complex_of_real ` {0<..<1 / a}"
        by (intro islimpt_isCont_image continuous_intros open_imp_islimpt)
           (use ab in auto simp: complex_eq_iff eventually_at_filter field_simps)
    next
      fix b assume "b  complex_of_real ` {0<..<1/a}"
      then obtain x where x: "b = complex_of_real x" "x  {0<..<1/a}"
        by auto
      show "f b = 0"
        unfolding f_def using real_eq1[of a x] x ab
        by (simp add: field_simps ramanujan_theta_of_real flip: qpochhammer_inf_of_real)
    qed (use ab in auto simp: norm_mult field_simps)
    thus ?thesis
      by (simp add: f_def)
  qed

  show "ramanujan_theta a b = (-a ; a * b) * (-b ; a * b) * (a * b ; a * b)"
  proof (cases "b = 0")
    case True
    thus ?thesis by auto
  next
    case False
    note ab = assms
    define f :: "complex  complex"
      where "f = (λa. ramanujan_theta a b - (-a ; a * b) * (-b ; a * b) * (a * b ; a * b))"

    have "f a = 0"
    proof (rule analytic_continuation[where f = f])
      show "f holomorphic_on (ball 0 (1 / norm b))"
        unfolding f_def using ab b  0
        by (intro holomorphic_intros) (auto simp: field_simps norm_mult)
    next
      show "complex_of_real ` {0<..<1/norm b}  ball 0 (1 / norm b)"
       and "complex_of_real (1/(2*norm b))  ball 0 (1 / norm b)"
        using ab b  0 by (auto simp: norm_mult norm_divide field_simps)
    next
      show "complex_of_real (1 / (2 * norm b)) islimpt
            complex_of_real ` {0<..<1 / norm b}"
        by (intro islimpt_isCont_image continuous_intros open_imp_islimpt)
           (use ab b  0 in auto simp: complex_eq_iff eventually_at_filter field_simps)
    next
      fix a assume "a  complex_of_real ` {0<..<1/norm b}"
      then obtain x where x: "a = complex_of_real x" "x  {0<..<1/norm b}"
        by auto
      show "f a = 0"
        unfolding f_def using real_eq2[of x b] x ab b  0
        by (simp add: norm_mult field_simps ramanujan_theta_of_real flip: qpochhammer_inf_of_real)
    qed (use ab b  0 in auto simp: norm_mult field_simps)
    thus ?thesis
      by (simp add: f_def)
  qed
qed

lemma ramanujan_theta_triple_product_real:
  fixes a b :: real
  assumes ab: "¦a * b¦ < 1"
  shows   "ramanujan_theta a b = (-a ; a * b) * (-b ; a * b) * (a * b ; a * b)"
proof -
  define a' b' where "a' = complex_of_real a" and "b' = complex_of_real b"
  have "complex_of_real (ramanujan_theta a b) = ramanujan_theta a' b'"
    by (simp add: a'_def b'_def ramanujan_theta_of_real)
  also have " = (-a' ; a' * b') * (-b' ; a' * b') * (a' * b' ; a' * b')"
    by (rule ramanujan_theta_triple_product_complex)
       (use assms in auto simp: a'_def b'_def abs_mult norm_mult)
  also have " = complex_of_real ((-a ; a * b) * (-b ; a * b) * (a * b ; a * b))"
    using assms by (simp flip: qpochhammer_inf_of_real add: a'_def b'_def)
  finally show ?thesis
    by (simp only: of_real_eq_iff)
qed 


subsection ‹The pentagonal number theorem›

text ‹
  An easy corollary of this is the Pentagonal Number Theorem, which, in our notation, simply reads:
  \[(q; q)_\infty = f(-q, -q^2) = \theta(-\sqrt{q}, q\sqrt{q})\]
  
›
corollary pentagonal_number_theorem_complex:
  fixes q :: complex
  assumes q: "norm q < 1"
  shows   "(q; q) = ramanujan_theta (-q) (-(q2))"
proof -
  have "ramanujan_theta (-q) (-(q2)) = (i<3. (q * q ^ i; q^3))"
    by (subst ramanujan_theta_triple_product_complex)
       (use q in simp_all flip: power_Suc add: norm_power power_less_one_iff eval_nat_numeral)
  also have " = (q ; q)"
    by (rule prod_qpochhammer_group) (use q in auto)
  finally show ?thesis ..
qed

lemma pentagonal_number_theorem_real:
  fixes q :: real
  assumes q: "¦q¦ < 1"
  shows   "(q; q) = ramanujan_theta (-q) (-(q2))"
proof -
  have "complex_of_real ((q; q)) = (of_real q; of_real q)"
    by (subst qpochhammer_inf_of_real) (use q in auto)
  also have " = complex_of_real (ramanujan_theta (-q) (-(q^2)))"
    by (subst pentagonal_number_theorem_complex)
       (use q in auto simp flip: ramanujan_theta_of_real)
  finally show ?thesis
    by (simp only: of_real_eq_iff)
qed

text ‹
  The following is the more explicit form of the Pentagonal Number Theorem usually found
  in textbooks:
  \[\prod_{n=1}^\infty \big(1 - q^n\big) = \sum_{k=-\infty}^\infty (-1)^k q^{k(3k-1)/2}\]
  The exponents $g_k = k(3k-1)/2$ (for $k\in\mathbb{Z}$) are called the \emph{generalised
  pentagonal numbers}.
›
corollary pentagonal_number_theorem_complex':
  fixes q :: complex
  assumes q: "norm q < 1"
  shows "abs_convergent_prod (λn. 1 - q^(n+1))" (is ?th1)
    and "(λk. (-1) powi k * q powi (k*(3*k-1) div 2)) abs_summable_on UNIV" (is ?th2)
    and "(n::nat. 1 - q^(n+1)) = ((k::int). (-1) powi k * q powi (k * (3*k-1) div 2))" (is ?th3)
proof -
  have "?th1  ?th2  ?th3"
  proof (cases "q = 0")
    case [simp]: True
    have "((λn. (-1) powi n * q powi (n*(3*n-1) div 2)) 
            has_sum ramanujan_theta (-q) (-q2)) {0}"
      by (intro has_sum_finiteI) auto
    also have "?this  ((λn. (-1) powi n * q powi (n*(3*n-1) div 2)) 
                         has_sum ramanujan_theta (-q) (-q2)) UNIV"
      by (intro has_sum_cong_neutral) (auto simp: dvd_div_eq_0_iff)
    finally show ?thesis
      by (auto simp: abs_convergent_prod_def has_sum_iff summable_on_iff_abs_summable_on_complex)
  next
    case [simp]: False
    have "(λn. 1 + norm q ^ Suc n) has_prod (-norm q; norm q)"
      using has_prod_qpochhammer_inf[of "norm q" "-norm q"] q by simp
    hence th1: "abs_convergent_prod (λn. 1 - q ^ (n+1))"
      by (simp add: abs_convergent_prod_def norm_mult norm_power has_prod_iff)
  
    have prod: "(λn. 1 - q ^ Suc n) has_prod (q; q)"
      using has_prod_qpochhammer_inf[of q q] q by simp
    have "((λn. (-q) powi (n*(n+1) div 2) * (-(q^2)) powi (n*(n-1) div 2))
             has_sum ramanujan_theta (-q) (-(q^2))) UNIV"
      by (rule has_sum_ramanujan_theta)
         (auto simp: norm_power power_less_one_iff q simp flip: power_Suc)
    also have "(λn. (-q) powi (n*(n+1) div 2) * (-(q^2)) powi (n*(n-1) div 2)) = 
               (λn. (- 1) powi n * q powi (n*(3*n-1) div 2))"
      (is "?lhs = ?rhs")
    proof
      fix n :: int
      have "(-q) powi (n*(n+1) div 2) * (-(q^2)) powi (n*(n-1) div 2) =
              (-1) powi (n*(n+1) div 2 + n*(n-1) div 2) * 
              (q powi (n*(n+1) div 2) * (q^2) powi (n*(n-1) div 2))"
        by (auto simp: power_int_minus_left)
      also have "n*(n+1) div 2 + n*(n-1) div 2 = (n*(n+1) + n*(n-1)) div 2"
        by (rule div_plus_div_distrib_dvd_left [symmetric]) auto
      also have "(n*(n+1) + n*(n-1)) div 2 = n ^ 2"
        by (simp add: algebra_simps power2_eq_square)
      also have "(-1) powi (n ^ 2) = (-1::complex) powi n"
        by (auto simp: power_int_minus_left)
      also have "(q^2) powi (n*(n-1) div 2) = q powi (n*(n-1))"
        by (simp add: power_int_power)
      also have "q powi (n * (n + 1) div 2) * q powi (n * (n - 1)) =
                 q powi (n * (n + 1) div 2 + (2 * n * (n - 1)) div 2)"
        by (simp add: power_int_add)
      also have "n * (n + 1) div 2 + (2 * n * (n - 1)) div 2 = (n*(n+1) + 2*n*(n-1)) div 2"
        by (rule div_plus_div_distrib_dvd_left [symmetric]) auto
      also have "n*(n+1) + 2*n*(n-1) = n * (3 * n - 1)"
        by (simp add: algebra_simps)
      finally show "?lhs n = ?rhs n" .
    qed
    finally have sum: "((λn. (-1) powi n * q powi (n*(3*n-1) div 2)) 
                         has_sum ramanujan_theta (-q) (-q2)) UNIV" .
  
    have "(λn. (-1) powi n * q powi (n*(3*n-1) div 2)) summable_on UNIV"
      using sum  by (rule has_sum_imp_summable)
    hence th2: "(λn. (-1) powi n * q powi (n*(3*n-1) div 2)) abs_summable_on UNIV"
      by (simp add: summable_on_iff_abs_summable_on_complex)

    have th3: "(n. 1 - q ^ (n+1)) = ((k::int). (-1) powi k * q powi (k * (3*k-1) div 2))"
      using sum prod pentagonal_number_theorem_complex[OF q]
      by (simp add: has_prod_iff has_sum_iff)
    
    show ?thesis
      using th1 th2 th3 by blast
  qed
  thus ?th1 ?th2 ?th3
    by blast+
qed


subsection ‹(Non-)vanishing of theta functinos›

text ‹
  A corollary of the Jacobi triple product: the Jacobi theta function has no zeros apart from
  the ``obvious'' ones, i.e.\ the ones at the center of the cells of the lattice generated
  by the period $1$ and the quasiperiod $t$.
›
corollary jacobi_theta_00_eq_0_iff_complex:
  fixes z t :: complex
  assumes "Im t > 0"
  shows   "jacobi_theta_00 z t = 0  (m n. z = (of_int m + 1/2) + (of_int n + 1/2) * t)"
proof
  assume "m n. z = (of_int m + 1/2) + (of_int n + 1/2) * t"
  then obtain m n where mn: "z = (of_int m + 1/2) + (of_int n + 1/2) * t"
    by blast
  show "jacobi_theta_00 z t = 0"
    unfolding mn by (rule jacobi_theta_00_eq_0')
next
  assume "jacobi_theta_00 z t = 0"
  define w q where "w = to_nome (2*z)" and "q = to_nome t"
  have [simp]: "w  0" "q  0"
    by (auto simp: w_def q_def)
  have q: "norm q < 1"
    using assms by (auto simp: q_def norm_to_nome)
  have q2: "norm (q ^ 2) < 1"
    by (simp add: norm_power power_less_one_iff q)

  note jacobi_theta_00 z t = 0
  also have "jacobi_theta_00 z t = jacobi_theta_nome w q"
    by (simp add: jacobi_theta_00_def w_def q_def to_nome_power)
  also have " = (q2 ; q2) * (- q * w ; q2) * (- q / w ; q2)"
    by (rule jacobi_theta_nome_triple_product_complex) (use q in auto)
  also have " = 0  ((n. q^(2*n+2) = 1)  (n. (w*q^(2*n+1)) = -1)  (n. (q^(2*n+1)/w) = -1))"
    using q q2 by (simp add: qpochhammer_inf_zero_iff power_add power_mult mult.assoc
                             power2_eq_square power_mult_distrib mult.commute[of w]
                             minus_equation_iff[of _ 1] eq_commute[of "-1"])
  also have "(λn. q^(2*n+2) = 1) = (λ_. False)"
  proof
    fix n :: nat
    have "norm (q ^ (2*n+2)) < 1"
      unfolding norm_power by (subst power_less_one_iff) (use q in auto)
    thus "q^(2*n+2) = 1  False"
      by auto
  qed
  finally show "m n. z = (of_int m + 1/2) + (of_int n + 1/2) * t"
  proof (elim disjE exE FalseE)
    fix n :: nat
    assume "w * q ^ (2*n+1) = -1"
    also have "w * q ^ (2*n+1) = to_nome (2 * z + (2 * of_nat n + 1) * t)"
      by (simp add: w_def q_def to_nome_power algebra_simps flip: to_nome_add)
    finally have "(2 * z + (2 * of_nat n + 1) * t - 1) / 2  "
      unfolding to_nome_eq_neg1_iff' .
    then obtain m where m: "(2 * z + (2 * of_nat n + 1) * t - 1) / 2 = of_int m"
      by (elim Ints_cases)
    have "z = (of_int m + 1/2) + (of_int (-int (n+1)) + 1/2) * t"
      by (subst m [symmetric]) (auto simp: field_simps)
    thus ?thesis
      by blast
  next
    fix n :: nat
    assume "q ^ (2*n+1) / w = -1"
    also have "q ^ (2*n+1) / w = to_nome (-2 * z + (2 * of_nat n + 1) * t)"
      by (simp add: w_def q_def to_nome_power algebra_simps flip: to_nome_add to_nome_diff)
    finally have "(-2 * z + (2 * of_nat n + 1) * t - 1) / 2  "
      unfolding to_nome_eq_neg1_iff' .
    then obtain m where m: "(-2 * z + (2 * of_nat n + 1) * t - 1) / 2 = of_int m"
      by (elim Ints_cases)
    have "z = (of_int (-m-1) + 1/2) + (of_int (int n) + 1/2) * t"
      unfolding of_int_diff of_int_minus by (subst m [symmetric]) (auto simp: field_simps)
    thus ?thesis
      by blast
  qed
qed

lemma jacobi_theta_00_nonzero:
  assumes z: "Im t > 0" and "Im z / Im t - 1 / 2  "
  shows   "jacobi_theta_00 z t  0"
proof
  assume "jacobi_theta_00 z t = 0"
  then obtain m n where *: "z = (of_int m + 1/2) + (of_int n + 1/2) * t"
    by (subst (asm) jacobi_theta_00_eq_0_iff_complex) (use Im t > 0 in auto)
  hence "Im z = (of_int n + 1/2) * Im t"
    by simp
  hence "Im z / Im t - 1 / 2 = of_int n"
    using z by (auto simp: field_simps)
  also have "  "
    by auto
  finally show False
    using assms by simp
qed

lemma jacobi_theta_00_0_left_nonzero:
  assumes "Im t > 0"
  shows   "jacobi_theta_00 0 t  0"
  by (rule jacobi_theta_00_nonzero) (use assms in auto)

lemma jacobi_theta_nome_nonzero_complex:
  fixes q w :: complex
  assumes [simp]: "w  0" "norm q < 1"
  assumes "q = 0  (ln (norm w) / ln (norm q) - 1) / 2  "
  shows   "jacobi_theta_nome w q  0"
proof (cases "q = 0")
  case [simp]: False
  define z where "z = -𝗂 * ln w / pi" 
  define t where "t = -𝗂 * ln q / pi" 
  have [simp]: "to_nome z = w" "to_nome t = q"
    using assms by (simp_all add: z_def t_def to_nome_def)
  have "Im t > 0"
    by (auto simp: t_def field_simps)

  have *: "Im z / (2 * Im t) - 1 / 2 = (ln (norm w) / ln (norm q) - 1) / 2"
    by (auto simp: z_def t_def)
  have "jacobi_theta_nome w q = jacobi_theta_00 (z/2) t"
    by (simp add: jacobi_theta_00_def to_nome_power)
  also have "  0"
  proof (rule jacobi_theta_00_nonzero)
    have "Im (z / 2) / Im t - 1/2 = ln (norm w) / (2 * ln (norm q)) - 1 / 2"
      by (simp add: z_def t_def)
    also have " = (ln (norm w) / ln (norm q) - 1) / 2"
      by (simp add: field_simps)
    also have "  "
      using assms by auto
    finally show "Im (z / 2) / Im t - 1/2  " .
  qed (use * Im t > 0 in auto)
  finally show ?thesis .
qed auto

lemma jacobi_theta_nome_q_q_nonzero_complex:
  assumes "norm (q::complex) < 1" "q  0"
  shows   "jacobi_theta_nome q q  0"
proof
  assume "jacobi_theta_nome q q = 0"
  define t where "t = -𝗂 * ln q / pi" 
  have [simp]: "to_nome t = q"
    using assms by (simp_all add: t_def to_nome_def)
  have "Im t > 0"
    using assms by (auto simp: t_def field_simps)

  note jacobi_theta_nome q q = 0
  also have "jacobi_theta_nome q q = jacobi_theta_00 (t/2) t"
    by (simp add: jacobi_theta_00_def to_nome_power)
  finally obtain m n where *: "t / 2 = (of_int m + 1/2) + (of_int n + 1/2) * t"
    by (subst (asm) jacobi_theta_00_eq_0_iff_complex) (use Im t > 0 in auto)
  have "t = (2 * of_int m + 1) + (2 * of_int n + 1) * t"
    using arg_cong[OF *, of "λx. x * 2"] by (simp add: ring_distribs mult_ac)
  hence **: "2 * of_int n * t = -(2 * of_int m + 1)"
    by (Groebner_Basis.algebra)
  have "n = 0"
    using Im t > 0 arg_cong[OF **, of Im] by simp
  with ** have "-2 * complex_of_int m = of_int 1"
    by simp
  also have "-2 * complex_of_int m = of_int (-2 * m)"
    by simp
  finally have "-2 * m = 1"
    by (simp only: of_int_eq_iff)
  thus False
    by presburger
qed

lemma jacobi_theta_nome_nonzero_real:
  fixes q w :: real
  assumes [simp]: "w  0" "norm q < 1" and "(ln ¦w¦ / ln ¦q¦ - 1) / 2  "
  shows   "jacobi_theta_nome w q  0"
proof -
  have "jacobi_theta_nome (complex_of_real w) (complex_of_real q)  0"
    by (rule jacobi_theta_nome_nonzero_complex) (use assms in auto)
  thus ?thesis
    by (simp add: jacobi_theta_nome_of_real)
qed

lemma jacobi_theta_nome_1_left_nonzero_complex:
  assumes "norm (q :: complex) < 1"
  shows   "jacobi_theta_nome 1 q  0"
  by (rule jacobi_theta_nome_nonzero_complex) (use assms in auto)

lemma jacobi_theta_nome_1_left_nonzero_real:
  assumes "¦q::real¦ < 1"
  shows   "jacobi_theta_nome 1 q  0"
  by (metis assms jacobi_theta_nome_1_left_nonzero_complex jacobi_theta_nome_of_real 
            norm_of_real of_real_0 of_real_1)

unbundle no_qpochhammer_inf_notation

end