(* File: Bernoulli_FPS.thy Author: Manuel Eberl <manuel@pruvisto.org> Connection of Bernoulli numbers to formal power series; proof B_n = 0 for odd n > 1; Akiyama-Tanigawa algorithm. *) section ‹Connection of Bernoulli numbers to formal power series› theory Bernoulli_FPS imports Bernoulli "HOL-Computational_Algebra.Computational_Algebra" "HOL-Combinatorics.Stirling" "HOL-Number_Theory.Number_Theory" begin subsection ‹Preliminaries› context factorial_semiring begin lemma multiplicity_prime_prime: "prime p ⟹ prime q ⟹ multiplicity p q = (if p = q then 1 else 0)" by (simp add: prime_multiplicity_other) lemma prime_prod_dvdI: fixes f :: "'b ⇒ 'a" assumes "finite A" assumes "⋀x. x ∈ A ⟹ prime (f x)" assumes "⋀x. x ∈ A ⟹ f x dvd y" assumes "inj_on f A" shows "prod f A dvd y" proof (cases "y = 0") case False have nz: "f x ≠ 0" if "x ∈ A" for x using assms(2)[of x] that by auto have "prod f A ≠ 0" using assms nz by (subst prod_zero_iff) auto thus ?thesis proof (rule multiplicity_le_imp_dvd) fix p :: 'a assume "prime p" show "multiplicity p (prod f A) ≤ multiplicity p y" proof (cases "p dvd prod f A") case True then obtain x where x: "x ∈ A" and "p dvd f x" using ‹prime p› assms by (subst (asm) prime_dvd_prod_iff) auto have "multiplicity p (prod f A) = (∑x∈A. multiplicity p (f x))" using assms ‹prime p› nz by (intro prime_elem_multiplicity_prod_distrib) auto also have "… = (∑x∈{x}. 1 :: nat)" using assms ‹prime p› ‹p dvd f x› primes_dvd_imp_eq x by (intro Groups_Big.sum.mono_neutral_cong_right) (auto simp: multiplicity_prime_prime inj_on_def) finally have "multiplicity p (prod f A) = 1" by simp also have "1 ≤ multiplicity p y" using assms nz ‹prime p› ‹y ≠ 0› x ‹p dvd f x› by (intro multiplicity_geI) force+ finally show ?thesis . qed (auto simp: not_dvd_imp_multiplicity_0) qed qed auto end (* TODO: Move? *) context semiring_gcd begin lemma gcd_add_dvd_right1: "a dvd b ⟹ gcd a (b + c) = gcd a c" by (elim dvdE) (simp add: gcd_add_mult mult.commute[of a]) lemma gcd_add_dvd_right2: "a dvd c ⟹ gcd a (b + c) = gcd a b" using gcd_add_dvd_right1[of a c b] by (simp add: add_ac) lemma gcd_add_dvd_left1: "a dvd b ⟹ gcd (b + c) a = gcd c a" using gcd_add_dvd_right1[of a b c] by (simp add: gcd.commute) lemma gcd_add_dvd_left2: "a dvd c ⟹ gcd (b + c) a = gcd b a" using gcd_add_dvd_right2[of a c b] by (simp add: gcd.commute) end context ring_gcd begin lemma gcd_diff_dvd_right1: "a dvd b ⟹ gcd a (b - c) = gcd a c" using gcd_add_dvd_right1[of a b "-c"] by simp lemma gcd_diff_dvd_right2: "a dvd c ⟹ gcd a (b - c) = gcd a b" using gcd_add_dvd_right2[of a "-c" b] by simp lemma gcd_diff_dvd_left1: "a dvd b ⟹ gcd (b - c) a = gcd c a" using gcd_add_dvd_left1[of a b "-c"] by simp lemma gcd_diff_dvd_left2: "a dvd c ⟹ gcd (b - c) a = gcd b a" using gcd_add_dvd_left2[of a "-c" b] by simp end lemma cong_int: "[a = b] (mod m) ⟹ [int a = int b] (mod m)" by (simp add: cong_int_iff) lemma Rats_int_div_natE: assumes "(x :: 'a :: field_char_0) ∈ ℚ" obtains m :: int and n :: nat where "n > 0" and "x = of_int m / of_nat n" and "coprime m n" proof - from assms obtain r where [simp]: "x = of_rat r" by (auto simp: Rats_def) obtain a b where [simp]: "r = Rat.Fract a b" and ab: "b > 0" "coprime a b" by (cases r) from ab show ?thesis by (intro that[of "nat b" a]) (auto simp: of_rat_rat) qed lemma sum_in_Ints: "(⋀x. x ∈ A ⟹ f x ∈ ℤ) ⟹ sum f A ∈ ℤ" by (induction A rule: infinite_finite_induct) auto lemma Ints_real_of_nat_divide: "b dvd a ⟹ real a / real b ∈ ℤ" by auto lemma product_dvd_fact: assumes "a > 1" "b > 1" "a = b ⟶ a > 2" shows "(a * b) dvd fact (a * b - 1)" proof (cases "a = b") case False have "a * 1 < a * b" and "1 * b < a * b" using assms by (intro mult_strict_left_mono mult_strict_right_mono; simp)+ hence ineqs: "a ≤ a * b - 1" "b ≤ a * b - 1" by linarith+ from False have "a * b = ∏{a,b}" by simp also have "… dvd ∏{1..a * b - 1}" using assms ineqs by (intro prod_dvd_prod_subset) auto finally show ?thesis by (simp add: fact_prod) next case [simp]: True from assms have "a > 2" by auto hence "a * 2 < a * b" using assms by (intro mult_strict_left_mono; simp) hence *: "2 * a ≤ a * b - 1" by linarith have "a * a dvd (2 * a) * a" by simp also have "… = ∏{2*a, a}" using assms by auto also have "… dvd ∏{1..a * b - 1}" using assms * by (intro prod_dvd_prod_subset) auto finally show ?thesis by (simp add: fact_prod) qed lemma composite_imp_factors_nat: assumes "m > 1" "¬prime (m::nat)" shows "∃n k. m = n * k ∧ 1 < n ∧ n < m ∧ 1 < k ∧ k < m" proof - from assms have "¬irreducible m" by (simp flip: prime_elem_iff_irreducible ) then obtain a where a: "a dvd m" "¬m dvd a" "a ≠ 1" using assms by (auto simp: irreducible_altdef) then obtain b where [simp]: "m = a * b" by auto from a assms have "a ≠ 0" "b ≠ 0" "b ≠ 1" by (auto intro!: Nat.gr0I) with a have "a > 1" "b > 1" by linarith+ moreover from this and a have "a < m" "b < m" by auto ultimately show ?thesis using ‹m = a * b› by blast qed text ‹ This lemma describes what the numerator and denominator of a finite subseries of the harmonic series are when it is written as a single fraction. › lemma sum_inverses_conv_fraction: fixes f :: "'a ⇒ 'b :: field" assumes "⋀x. x ∈ A ⟹ f x ≠ 0" "finite A" shows "(∑x∈A. 1 / f x) = (∑x∈A. ∏y∈A-{x}. f y) / (∏x∈A. f x)" proof - have "(∑x∈A. (∏y∈A. f y) / f x) = (∑x∈A. ∏y∈A-{x}. f y)" using prod.remove[of A _ f] assms by (intro sum.cong refl) (auto simp: field_simps) thus ?thesis using assms by (simp add: field_simps sum_distrib_right sum_distrib_left) qed text ‹ If all terms in the subseries are primes, this fraction is automatically on lowest terms. › lemma sum_prime_inverses_fraction_coprime: fixes f :: "'a ⇒ nat" assumes "finite A" and primes: "⋀x. x ∈ A ⟹ prime (f x)" and inj: "inj_on f A" defines "a ≡ (∑x∈A. ∏y∈A-{x}. f y)" shows "coprime a (∏x∈A. f x)" proof (intro prod_coprime_right) fix x assume x: "x ∈ A" have "a = (∏y∈A-{x}. f y) + (∑y∈A-{x}. ∏z∈A-{y}. f z)" unfolding a_def using ‹finite A› and x by (rule sum.remove) also have "gcd … (f x) = gcd (∏y∈A-{x}. f y) (f x)" using ‹finite A› and x by (intro gcd_add_dvd_left2 dvd_sum dvd_prodI) auto also from x primes inj have "coprime (∏y∈A-{x}. f y) (f x)" by (intro prod_coprime_left) (auto intro!: primes_coprime simp: inj_on_def) hence "gcd (∏y∈A-{x}. f y) (f x) = 1" by simp finally show "coprime a (f x)" by (simp only: coprime_iff_gcd_eq_1) qed (* END TODO *) text ‹ In the following, we will prove the correctness of the Akiyama--Tanigawa algorithm~\<^cite>‹"kaneko2000"›, which is a simple algorithm for computing Bernoulli numbers that was discovered by Akiyama and Tanigawa~\<^cite>‹"aki_tani1999"› essentially as a by-product of their studies of the Euler--Zagier multiple zeta function. The algorithm is based on a number triangle (similar to Pascal's triangle) in which the Bernoulli numbers are the leftmost diagonal. While the algorithm itself is quite simple, proving it correct is not entirely trivial. We will use generating functions and Stirling numbers, mostly following the presentation by Kaneko~\<^cite>‹"kaneko2000"›. › text ‹ The following operator is a variant of the @{term fps_XD} operator where the multiplication is not with @{term fps_X}, but with an arbitrary formal power series. It is not quite clear if this operator has a less ad-hoc meaning than the fashion in which we use it; it is, however, very useful for proving the relationship between Stirling numbers and Bernoulli numbers. › context includes fps_notation begin definition fps_XD' where "fps_XD' a = (λb. a * fps_deriv b)" lemma fps_XD'_0 [simp]: "fps_XD' a 0 = 0" by (simp add: fps_XD'_def) lemma fps_XD'_1 [simp]: "fps_XD' a 1 = 0" by (simp add: fps_XD'_def) lemma fps_XD'_fps_const [simp]: "fps_XD' a (fps_const b) = 0" by (simp add: fps_XD'_def) lemma fps_XD'_fps_of_nat [simp]: "fps_XD' a (of_nat b) = 0" by (simp add: fps_XD'_def) lemma fps_XD'_fps_of_int [simp]: "fps_XD' a (of_int b) = 0" by (simp add: fps_XD'_def) lemma fps_XD'_fps_numeral [simp]: "fps_XD' a (numeral b) = 0" by (simp add: fps_XD'_def) lemma fps_XD'_add [simp]: "fps_XD' a (b + c :: 'a :: comm_ring_1 fps) = fps_XD' a b + fps_XD' a c" by (simp add: fps_XD'_def algebra_simps) lemma fps_XD'_minus [simp]: "fps_XD' a (b - c :: 'a :: comm_ring_1 fps) = fps_XD' a b - fps_XD' a c" by (simp add: fps_XD'_def algebra_simps) lemma fps_XD'_prod: "fps_XD' a (b * c :: 'a :: comm_ring_1 fps) = fps_XD' a b * c + b * fps_XD' a c" by (simp add: fps_XD'_def algebra_simps) lemma fps_XD'_power: "fps_XD' a (b ^ n :: 'a :: idom fps) = of_nat n * b ^ (n - 1) * fps_XD' a b" proof (cases "n = 0") case False have "b * fps_XD' a (b ^ n) = of_nat n * b ^ n * fps_XD' a b" by (induction n) (simp_all add: fps_XD'_prod algebra_simps) also have "… = b * (of_nat n * b ^ (n - 1) * fps_XD' a b)" by (cases n) (simp_all add: algebra_simps) finally show ?thesis using False by (subst (asm) mult_cancel_left) (auto simp: power_0_left) qed simp_all lemma fps_XD'_power_Suc: "fps_XD' a (b ^ Suc n :: 'a :: idom fps) = of_nat (Suc n) * b ^ n * fps_XD' a b" by (subst fps_XD'_power) simp_all lemma fps_XD'_sum: "fps_XD' a (sum f A) = sum (λx. fps_XD' (a :: 'a :: comm_ring_1 fps) (f x)) A" by (induction A rule: infinite_finite_induct) simp_all lemma fps_XD'_funpow_affine: fixes G H :: "real fps" assumes [simp]: "fps_deriv G = 1" defines "S ≡ λn i. fps_const (real (Stirling n i))" shows "(fps_XD' G ^^ n) H = (∑m≤n. S n m * G ^ m * (fps_deriv ^^ m) H)" proof (induction n arbitrary: H) case 0 thus ?case by (simp add: S_def) next case (Suc n H) have "(∑m≤Suc n. S (Suc n) m * G ^ m * (fps_deriv ^^ m) H) = (∑i≤n. of_nat (Suc i) * S n (Suc i) * G ^ Suc i * (fps_deriv ^^ Suc i) H) + (∑i≤n. S n i * G ^ Suc i * (fps_deriv ^^ Suc i) H)" (is "_ = sum (λi. ?f (Suc i)) … + ?S2") by (subst sum.atMost_Suc_shift) (simp_all add: sum.distrib algebra_simps fps_of_nat S_def fps_const_add [symmetric] fps_const_mult [symmetric] del: fps_const_add fps_const_mult) also have "sum (λi. ?f (Suc i)) {..n} = sum (λi. ?f (Suc i)) {..<n}" by (intro sum.mono_neutral_right) (auto simp: S_def) also have "… = ?f 0 + …" by simp also have "… = sum ?f {..n}" by (subst sum.atMost_shift [symmetric]) simp_all also have "… + ?S2 = (∑x≤n. fps_XD' G (S n x * G ^ x * (fps_deriv ^^ x) H))" unfolding sum.distrib [symmetric] proof (rule sum.cong, goal_cases) case (2 i) thus ?case unfolding fps_XD'_prod fps_XD'_power by (cases i) (auto simp: fps_XD'_prod fps_XD'_power_Suc algebra_simps of_nat_diff S_def fps_XD'_def) qed simp_all also have "… = (fps_XD' G ^^ Suc n) H" by (simp add: Suc.IH fps_XD'_sum) finally show ?case .. qed subsection ‹Generating function of Stirling numbers› lemma Stirling_n_0: "Stirling n 0 = (if n = 0 then 1 else 0)" by (cases n) simp_all text ‹ The generating function of Stirling numbers w.\,r.\,t.\ their first argument: \[\sum_{n=0}^\infty \genfrac{\{}{\}}{0pt}{}{n}{m} \frac{x^n}{n!} = \frac{(e^x - 1)^m}{m!}\] › definition Stirling_fps :: "nat ⇒ real fps" where "Stirling_fps m = fps_const (1 / fact m) * (fps_exp 1 - 1) ^ m" theorem sum_Stirling_binomial: "Stirling (Suc n) (Suc m) = (∑i = 0..n. Stirling i m * (n choose i))" proof - have "real (Stirling (Suc n) (Suc m)) = real (∑i = 0..n. Stirling i m * (n choose i))" proof (induction n arbitrary: m) case (Suc n m) have "real (∑i = 0..Suc n. Stirling i m * (Suc n choose i)) = real (∑i = 0..n. Stirling (Suc i) m * (Suc n choose Suc i)) + real (Stirling 0 m)" by (subst sum.atLeast0_atMost_Suc_shift) simp_all also have "real (∑i = 0..n. Stirling (Suc i) m * (Suc n choose Suc i)) = real (∑i = 0..n. (n choose i) * Stirling (Suc i) m) + real (∑i = 0..n. (n choose Suc i) * Stirling (Suc i) m)" by (simp add: algebra_simps sum.distrib) also have "(∑i = 0..n. (n choose Suc i) * Stirling (Suc i) m) = (∑i = Suc 0..Suc n. (n choose i) * Stirling i m)" by (subst sum.shift_bounds_cl_Suc_ivl) simp_all also have "… = (∑i = Suc 0..n. (n choose i) * Stirling i m)" by (intro sum.mono_neutral_right) auto also have "… = real (∑i = 0..n. Stirling i m * (n choose i)) - real (Stirling 0 m)" by (simp add: sum.atLeast_Suc_atMost mult_ac) also have "real (∑i = 0..n. Stirling i m * (n choose i)) = real (Stirling (Suc n) (Suc m))" by (rule Suc.IH [symmetric]) also have "real (∑i = 0..n. (n choose i) * Stirling (Suc i) m) = real m * real (Stirling (Suc n) (Suc m)) + real (Stirling (Suc n) m)" by (cases m; (simp only: Suc.IH, simp add: algebra_simps sum.distrib sum_distrib_left sum_distrib_right)) also have "… + (real (Stirling (Suc n) (Suc m)) - real (Stirling 0 m)) + real (Stirling 0 m) = real (Suc m * Stirling (Suc n) (Suc m) + Stirling (Suc n) m)" by (simp add: algebra_simps del: Stirling.simps) also have "Suc m * Stirling (Suc n) (Suc m) + Stirling (Suc n) m = Stirling (Suc (Suc n)) (Suc m)" by (rule Stirling.simps(4) [symmetric]) finally show ?case .. qed simp_all thus ?thesis by (subst (asm) of_nat_eq_iff) qed lemma Stirling_fps_aux: "(fps_exp 1 - 1) ^ m $ n * fact n = fact m * real (Stirling n m)" proof (induction m arbitrary: n) case 0 thus ?case by (simp add: Stirling_n_0) next case (Suc m n) show ?case proof (cases n) case 0 thus ?thesis by simp next case (Suc n') hence "(fps_exp 1 - 1 :: real fps) ^ Suc m $ n * fact n = fps_deriv ((fps_exp 1 - 1) ^ Suc m) $ n' * fact n'" by (simp_all add: algebra_simps del: power_Suc) also have "fps_deriv ((fps_exp 1 - 1 :: real fps) ^ Suc m) = fps_const (real (Suc m)) * ((fps_exp 1 - 1) ^ m * fps_exp 1)" by (subst fps_deriv_power) simp_all also have "… $ n' * fact n' = real (Suc m) * ((∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i / fact (n' - i)) * fact n')" unfolding fps_mult_left_const_nth by (simp add: fps_mult_nth Suc.IH sum_distrib_right del: of_nat_Suc) also have "(∑i = 0..n'. (fps_exp 1 - 1 :: real fps) ^ m $ i / fact (n' - i)) * fact n' = (∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i * fact n' / fact (n' - i))" by (subst sum_distrib_right, rule sum.cong) (simp_all add: divide_simps) also have "… = (∑i = 0..n'. (fps_exp 1 - 1) ^ m $ i * fact i * (n' choose i))" by (intro sum.cong refl) (simp_all add: binomial_fact) also have "… = (∑i = 0..n'. fact m * real (Stirling i m) * real (n' choose i))" by (simp only: Suc.IH) also have "real (Suc m) * … = fact (Suc m) * (∑i = 0..n'. real (Stirling i m) * real (n' choose i))" (is "_ = _ * ?S") by (simp add: sum_distrib_left sum_distrib_right mult_ac del: of_nat_Suc) also have "?S = Stirling (Suc n') (Suc m)" by (subst sum_Stirling_binomial) simp also have "Suc n' = n" by (simp add: Suc) finally show ?thesis . qed qed lemma Stirling_fps_nth: "Stirling_fps m $ n = Stirling n m / fact n" unfolding Stirling_fps_def using Stirling_fps_aux[of m n] by (simp add: field_simps) theorem Stirling_fps_altdef: "Stirling_fps m = Abs_fps (λn. Stirling n m / fact n)" by (simp add: fps_eq_iff Stirling_fps_nth) theorem Stirling_closed_form: "real (Stirling n k) = (∑j≤k. (-1)^(k - j) * real (k choose j) * real j ^ n) / fact k" proof - have "(fps_exp 1 - 1 :: real fps) = (fps_exp 1 + (-1))" by simp also have "… ^ k = (∑j≤k. of_nat (k choose j) * fps_exp 1 ^ j * (- 1) ^ (k - j))" unfolding binomial_ring .. also have "… = (∑j≤k. fps_const ((-1) ^ (k - j) * real (k choose j)) * fps_exp (real j))" by (simp add: fps_const_mult [symmetric] fps_const_power [symmetric] fps_const_neg [symmetric] mult_ac fps_of_nat fps_exp_power_mult del: fps_const_mult fps_const_power fps_const_neg) also have "… $ n = (∑j≤k. (- 1) ^ (k - j) * real (k choose j) * real j ^ n) / fact n" by (simp add: fps_sum_nth sum_divide_distrib) also have "… * fact n = (∑j≤k. (- 1) ^ (k - j) * real (k choose j) * real j ^ n)" by simp also note Stirling_fps_aux[of k n] finally show ?thesis by (simp add: atLeast0AtMost field_simps) qed subsection ‹Generating function of Bernoulli numbers› text ‹ We will show that the negative and positive Bernoulli numbers are the coefficients of the exponential generating function $\frac{x}{e^x - 1}$ (resp. $\frac{x}{1-e^{-x}}$), i.\,e. \[\sum_{n=0}^\infty B_n^{-} \frac{x^n}{n!} = \frac{x}{e^x - 1}\] \[\sum_{n=0}^\infty B_n^{+} \frac{x^n}{n!} = \frac{x}{1 - e^{-1}}\] › definition bernoulli_fps :: "'a :: real_normed_field fps" where "bernoulli_fps = fps_X / (fps_exp 1 - 1)" definition bernoulli'_fps :: "'a :: real_normed_field fps" where "bernoulli'_fps = fps_X / (1 - (fps_exp (-1)))" lemma bernoulli_fps_altdef: "bernoulli_fps = Abs_fps (λn. of_real (bernoulli n) / fact n :: 'a)" and bernoulli_fps_aux: "bernoulli_fps * (fps_exp 1 - 1 :: 'a :: real_normed_field fps) = fps_X" proof - have *: "Abs_fps (λn. of_real (bernoulli n) / fact n :: 'a) * (fps_exp 1 - 1) = fps_X" proof (rule fps_ext) fix n have "(Abs_fps (λn. of_real (bernoulli n) / fact n :: 'a) * (fps_exp 1 - 1)) $ n = (∑i = 0..n. of_real (bernoulli i) * (1 / fact (n - i) - (if n = i then 1 else 0)) / fact i)" by (auto simp: fps_mult_nth divide_simps split: if_splits intro!: sum.cong) also have "… = (∑i = 0..n. of_real (bernoulli i) / (fact i * fact (n - i)) - (if n = i then of_real (bernoulli i) / fact i else 0))" by (intro sum.cong) (simp_all add: field_simps) also have "… = (∑i = 0..n. of_real (bernoulli i) / (fact i * fact (n - i))) - of_real (bernoulli n) / fact n" unfolding sum_subtractf by (subst sum.delta') simp_all also have "… = (∑i<n. of_real (bernoulli i) / (fact i * fact (n - i)))" by (cases n) (simp_all add: atLeast0AtMost lessThan_Suc_atMost [symmetric]) also have "… = (∑i<n. fact n * (of_real (bernoulli i) / (fact i * fact (n - i)))) / fact n" by (subst sum_distrib_left [symmetric]) simp_all also have "(∑i<n. fact n * (of_real (bernoulli i) / (fact i * fact (n - i)))) = (∑i<n. of_nat (n choose i) * of_real (bernoulli i) :: 'a)" by (intro sum.cong) (simp_all add: binomial_fact) also have "… = of_real (∑i<n. (n choose i) * bernoulli i)" by simp also have "… / fact n = fps_X $ n" by (subst sum_binomial_times_bernoulli') simp_all finally show "(Abs_fps (λn. of_real (bernoulli n) / fact n :: 'a) * (fps_exp 1 - 1)) $ n = fps_X $ n" . qed moreover show "bernoulli_fps = Abs_fps (λn. of_real (bernoulli n) / fact n :: 'a)" unfolding bernoulli_fps_def by (subst * [symmetric]) simp_all ultimately show "bernoulli_fps * (fps_exp 1 - 1 :: 'a fps) = fps_X" by simp qed theorem fps_nth_bernoulli_fps [simp]: "fps_nth bernoulli_fps n = of_real (bernoulli n) / fact n" by (simp add: bernoulli_fps_altdef) lemma bernoulli'_fps_aux: "(fps_exp 1 - 1) * Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a) = fps_exp 1 * fps_X" and bernoulli'_fps_aux': "(1 - fps_exp (-1)) * Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a) = fps_X" and bernoulli'_fps_altdef: "bernoulli'_fps = Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a :: real_normed_field)" proof - have "Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a) = bernoulli_fps + fps_X" by (simp add: fps_eq_iff bernoulli'_def) also have "(fps_exp 1 - 1) * … = fps_exp 1 * fps_X" using bernoulli_fps_aux by (simp add: algebra_simps) finally show "(fps_exp 1 - 1) * Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a) = fps_exp 1 * fps_X" . also have "(fps_exp 1 - 1) = fps_exp 1 * (1 - fps_exp (-1 :: 'a))" by (simp add: algebra_simps fps_exp_add_mult [symmetric]) also note mult.assoc finally show *: "(1 - fps_exp (-1)) * Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a) = fps_X" by (subst (asm) mult_left_cancel) simp_all show "bernoulli'_fps = Abs_fps (λn. of_real (bernoulli' n) / fact n :: 'a)" unfolding bernoulli'_fps_def by (subst * [symmetric]) simp_all qed theorem fps_nth_bernoulli'_fps [simp]: "fps_nth bernoulli'_fps n = of_real (bernoulli' n) / fact n" by (simp add: bernoulli'_fps_altdef) lemma bernoulli_fps_conv_bernoulli'_fps: "bernoulli_fps = bernoulli'_fps - fps_X" by (simp add: fps_eq_iff bernoulli'_def) lemma bernoulli'_fps_conv_bernoulli_fps: "bernoulli'_fps = bernoulli_fps + fps_X" by (simp add: fps_eq_iff bernoulli'_def) theorem bernoulli_odd_eq_0: assumes "n ≠ 1" and "odd n" shows "bernoulli n = 0" proof - from bernoulli_fps_aux have "2 * bernoulli_fps * (fps_exp 1 - 1) = 2 * fps_X" by simp hence "(2 * bernoulli_fps + fps_X) * (fps_exp 1 - 1) = fps_X * (fps_exp 1 + 1)" by (simp add: algebra_simps) also have "fps_exp 1 - 1 = fps_exp (1/2) * (fps_exp (1/2) - fps_exp (-1/2 :: real))" by (simp add: algebra_simps fps_exp_add_mult [symmetric]) also have "fps_exp 1 + 1 = fps_exp (1/2) * (fps_exp (1/2) + fps_exp (-1/2 :: real))" by (simp add: algebra_simps fps_exp_add_mult [symmetric]) finally have "fps_exp (1/2) * ((2 * bernoulli_fps + fps_X) * (fps_exp (1/2) - fps_exp (- 1/2))) = fps_exp (1/2) * (fps_X * (fps_exp (1/2) + fps_exp (-1/2 :: real)))" by (simp add: algebra_simps) hence *: "(2 * bernoulli_fps + fps_X) * (fps_exp (1/2) - fps_exp (- 1/2)) = fps_X * (fps_exp (1/2) + fps_exp (-1/2 :: real))" (is "?lhs = ?rhs") by (subst (asm) mult_cancel_left) simp_all have "fps_compose ?lhs (-fps_X) = fps_compose ?rhs (-fps_X)" by (simp only: *) also have "fps_compose ?lhs (-fps_X) = (-2 * (bernoulli_fps oo - fps_X) + fps_X) * (fps_exp ((1/2)) - fps_exp (-1/2))" by (simp add: fps_compose_mult_distrib fps_compose_add_distrib fps_compose_sub_distrib algebra_simps) also have "fps_compose ?rhs (-fps_X) = -?rhs" by (simp add: fps_compose_mult_distrib fps_compose_add_distrib fps_compose_sub_distrib) also note * [symmetric] also have "- ((2 * bernoulli_fps + fps_X) * (fps_exp (1/2) - fps_exp (-1/2))) = ((-2 * bernoulli_fps - fps_X) * (fps_exp (1/2) - fps_exp (-1/2)))" by (simp add: algebra_simps) finally have "2 * (bernoulli_fps oo - fps_X) = 2 * (bernoulli_fps + fps_X :: real fps)" by (subst (asm) mult_cancel_right) (simp add: algebra_simps) hence **: "bernoulli_fps oo -fps_X = (bernoulli_fps + fps_X :: real fps)" by (subst (asm) mult_cancel_left) simp from assms have "(bernoulli_fps oo -fps_X) $ n = bernoulli n / fact n" by (subst **) simp also have "-fps_X = fps_const (-1 :: real) * fps_X" by (simp only: fps_const_neg [symmetric] fps_const_1_eq_1) simp also from assms have "(bernoulli_fps oo …) $ n = - bernoulli n / fact n" by (subst fps_compose_linear) simp finally show ?thesis by simp qed lemma bernoulli'_odd_eq_0: "n ≠ 1 ⟹ odd n ⟹ bernoulli' n = 0" by (simp add: bernoulli'_def bernoulli_odd_eq_0) text ‹ The following simplification rule takes care of rewriting @{term "bernoulli n"} to $0$ for any odd numeric constant greater than $1$: › lemma bernoulli_odd_numeral_eq_0 [simp]: "bernoulli (numeral (Num.Bit1 n)) = 0" by (rule bernoulli_odd_eq_0[OF _ odd_numeral]) auto lemma bernoulli'_odd_numeral_eq_0 [simp]: "bernoulli' (numeral (Num.Bit1 n)) = 0" by (simp add: bernoulli'_def) text ‹ The following explicit formula for Bernoulli numbers can also derived reasonably easily using the generating functions of Stirling numbers and Bernoulli numbers. The proof follows an answer by Marko Riedel on the Mathematics StackExchange~\<^cite>‹"riedel_mathse_2014"›. › theorem bernoulli_altdef: "bernoulli n = (∑m≤n. ∑k≤m. (-1)^k * real (m choose k) * real k^n / real (Suc m))" proof - have "(∑m≤n. ∑k≤m. (-1)^k * real (m choose k) * real k^n / real (Suc m)) = (∑m≤n. (∑k≤m. (-1)^k * real (m choose k) * real k^n) / real (Suc m))" by (subst sum_divide_distrib) simp_all also have "… = fact n * (∑m≤n. (- 1) ^ m / real (Suc m) * (fps_exp 1 - 1) ^ m $ n)" proof (subst sum_distrib_left, intro sum.cong refl) fix m assume m: "m ∈ {..n}" have "(∑k≤m. (-1)^k * real (m choose k) * real k^n) = (-1)^m * (∑k≤m. (-1)^(m - k) * real (m choose k) * real k^n)" by (subst sum_distrib_left, intro sum.cong refl) (auto simp: minus_one_power_iff) also have "… = (-1) ^ m * (real (Stirling n m) * fact m)" by (subst Stirling_closed_form) simp_all also have "real (Stirling n m) = Stirling_fps m $ n * fact n" by (subst Stirling_fps_nth) simp_all also have "… * fact m = (fps_exp 1 - 1) ^ m $ n * fact n" by (simp add: Stirling_fps_def) finally show "(∑k≤m. (-1)^k * real (m choose k) * real k^n) / real (Suc m) = fact n * ((- 1) ^ m / real (Suc m) * (fps_exp 1 - 1) ^ m $ n)" by simp qed also have "(∑m≤n. (- 1) ^ m / real (Suc m) * (fps_exp 1 - 1) ^ m $ n) = fps_compose (Abs_fps (λm. (-1) ^ m / real (Suc m))) (fps_exp 1 - 1) $ n" by (simp add: fps_compose_def atLeast0AtMost fps_sum_nth) also have "fps_ln 1 = fps_X * Abs_fps (λm. (-1) ^ m / real (Suc m))" unfolding fps_ln_def by (auto simp: fps_eq_iff) hence "Abs_fps (λm. (-1) ^ m / real (Suc m)) = fps_ln 1 / fps_X" by (metis fps_X_neq_zero nonzero_mult_div_cancel_left) also have "fps_compose … (fps_exp 1 - 1) = fps_compose (fps_ln 1) (fps_exp 1 - 1) / (fps_exp 1 - 1)" by (subst fps_compose_divide_distrib) auto also have "fps_compose (fps_ln 1) (fps_exp 1 - 1 :: real fps) = fps_X" by (simp add: fps_ln_fps_exp_inv fps_inv_fps_exp_compose) also have "(fps_X / (fps_exp 1 - 1)) = bernoulli_fps" by (simp add: bernoulli_fps_def) also have "fact n * … $ n = bernoulli n" by simp finally show ?thesis .. qed corollary%important bernoulli_conv_Stirling: "bernoulli n = (∑k≤n. (-1) ^ k * fact k / real (k + 1) * Stirling n k)" proof - have "(∑k≤n. (-1) ^ k * fact k / (k + 1) * Stirling n k) = (∑k≤n. ∑i≤k. (-1) ^ i * (k choose i) * i ^ n / real (k + 1))" proof (intro sum.cong, goal_cases) case (2 k) have "(-1) ^ k * fact k / (k + 1) * Stirling n k = (∑j≤k. (-1) ^ k * (-1) ^ (k - j) * (k choose j) * j ^ n / (k + 1))" by (simp add: Stirling_closed_form sum_distrib_left sum_divide_distrib mult_ac) also have "… = (∑j≤k. (-1) ^ j * (k choose j) * j ^ n / (k + 1))" by (intro sum.cong) (auto simp: uminus_power_if split: if_splits) finally show ?case . qed auto also have "… = bernoulli n" by (simp add: bernoulli_altdef) finally show ?thesis .. qed subsection ‹Von Staudt--Clausen Theorem› lemma vonStaudt_Clausen_lemma: assumes "n > 0" and "prime p" shows "[(∑m<p. (-1) ^ m * ((p - 1) choose m) * m ^ (2*n)) = (if (p - 1) dvd (2 * n) then -1 else 0)] (mod p)" proof (cases "(p - 1) dvd (2 * n)") case True have cong_power_2n: "[m ^ (2 * n) = 1] (mod p)" if "m > 0" "m < p" for m proof - from True obtain q where "2 * n = (p - 1) * q" by blast hence "[m ^ (2 * n) = (m ^ (p - 1)) ^ q] (mod p)" by (simp add: power_mult) also have "[(m ^ (p - 1)) ^ q = 1 ^ q] (mod p)" using assms ‹m > 0› ‹m < p› by (intro cong_pow fermat_theorem) auto finally show ?thesis by simp qed have "(∑m<p. (-1)^m * ((p - 1) choose m) * m ^ (2*n)) = (∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * m ^ (2*n))" using ‹n > 0› by (intro sum.mono_neutral_right) auto also have "[… = (∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * int 1)] (mod p)" by (intro cong_sum cong_mult cong_power_2n cong_int) auto also have "(∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * int 1) = (∑m∈insert 0 {0<..<p}. (-1)^m * ((p - 1) choose m)) - 1" by (subst sum.insert) auto also have "insert 0 {0<..<p} = {..p-1}" using assms prime_gt_0_nat[of p] by auto also have "(∑m≤p-1. (-1)^m * ((p - 1) choose m)) = 0" using prime_gt_1_nat[of p] assms by (subst choose_alternating_sum) auto finally show ?thesis using True by simp next case False define n' where "n' = (2 * n) mod (p - 1)" from assms False have "n' > 0" by (auto simp: n'_def dvd_eq_mod_eq_0) from False have "p ≠ 2" by auto with assms have "odd p" using prime_prime_factor two_is_prime_nat by blast have cong_pow_2n: "[m ^ (2*n) = m ^ n'] (mod p)" if "m > 0" "m < p" for m proof - from assms and that have "coprime p m" by (intro prime_imp_coprime) auto have "[2 * n = n'] (mod (p - 1))" by (simp add: n'_def) moreover have "ord p m dvd (p - 1)" using order_divides_totient[of p m] ‹coprime p m› assms by (auto simp: totient_prime) ultimately have "[2 * n = n'] (mod ord p m)" by (rule cong_dvd_modulus_nat) thus ?thesis using ‹coprime p m› by (subst order_divides_expdiff) auto qed have "(∑m<p. (-1)^m * ((p - 1) choose m) * m ^ (2*n)) = (∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * m ^ (2*n))" using ‹n > 0› by (intro sum.mono_neutral_right) auto also have "[… = (∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * m ^ n')] (mod p)" by (intro cong_sum cong_mult cong_pow_2n cong_int) auto also have "(∑m∈{0<..<p}. (-1)^m * ((p - 1) choose m) * m ^ n') = (∑m≤p-1. (-1)^m * ((p - 1) choose m) * m ^ n')" using ‹n' > 0› by (intro sum.mono_neutral_left) auto also have "… = (∑m≤p-1. (-1)^(p - Suc m) * ((p - 1) choose m) * m ^ n')" using ‹n' > 0› assms ‹odd p› by (intro sum.cong) (auto simp: uminus_power_if) also have "… = 0" proof - have "of_int (∑m≤p-1. (-1)^(p - Suc m) * ((p - 1) choose m) * m ^ n') = real (Stirling n' (p - 1)) * fact (p - 1)" by (simp add: Stirling_closed_form) also have "n' < p - 1" using assms prime_gt_1_nat[of p] by (auto simp: n'_def) hence "Stirling n' (p - 1) = 0" by simp finally show ?thesis by linarith qed finally show ?thesis using False by simp qed text ‹ The Von Staudt--Clausen theorem states that for ‹n > 0›, \[B_{2n} + \sum\limits_{p - 1\mid 2n} \frac{1}{p}\] is an integer. › theorem vonStaudt_Clausen: assumes "n > 0" shows "bernoulli (2 * n) + (∑p | prime p ∧ (p - 1) dvd (2 * n). 1 / real p) ∈ ℤ" (is "_ + ?P ∈ ℤ") proof - define P :: "nat ⇒ real" where "P = (λm. if prime (m + 1) ∧ m dvd (2 * n) then 1 / (m + 1) else 0)" define P' :: "nat ⇒ int" where "P' = (λm. if prime (m + 1) ∧ m dvd (2 * n) then 1 else 0)" have "?P = (∑p | prime (p + 1) ∧ p dvd (2 * n). 1 / real (p + 1))" by (rule sum.reindex_bij_witness[of _ "λp. p + 1" "λp. p - 1"]) (use prime_gt_0_nat in auto) also have "… = (∑m≤2*n. P m)" using ‹n > 0› by (intro sum.mono_neutral_cong_left) (auto simp: P_def dest!: dvd_imp_le) finally have "bernoulli (2 * n) + ?P = (∑m≤2*n. (-1)^m * (of_int (fact m * Stirling (2*n) m) / (m + 1)) + P m)" by (simp add: sum.distrib bernoulli_conv_Stirling sum_divide_distrib algebra_simps) also have "… = (∑m≤2*n. of_int ((-1)^m * fact m * Stirling (2*n) m + P' m) / (m + 1))" by (intro sum.cong) (auto simp: P'_def P_def field_simps) also have "… ∈ ℤ" proof (rule sum_in_Ints, goal_cases) case (1 m) have "m = 0 ∨ m = 3 ∨ prime (m + 1) ∨ (¬prime (m + 1) ∧ m > 3)" by (cases "m = 1"; cases "m = 2") (auto simp flip: numeral_2_eq_2) then consider "m = 0" | "m = 3" | "prime (m + 1)" | "¬prime (m + 1)" "m > 3" by blast thus ?case proof cases assume "m = 0" thus ?case by auto next assume [simp]: "m = 3" have "real_of_int (fact m * Stirling (2 * n) m) = real_of_int (9 ^ n + 3 - 3 * 4 ^ n)" using ‹n > 0› by (auto simp: P'_def fact_numeral Stirling_closed_form power_mult atMost_nat_numeral binomial_fact zero_power) hence "int (fact m * Stirling (2 * n) m) = 9 ^ n + 3 - 3 * 4 ^ n" by linarith also have "[… = 1 ^ n + (-1) - 3 * 0 ^ n] (mod 4)" by (intro cong_add cong_diff cong_mult cong_pow) (auto simp: cong_def) finally have dvd: "4 dvd int (fact m * Stirling (2 * n) m)" using ‹n > 0› by (simp add: cong_0_iff zero_power) have "real_of_int ((- 1) ^ m * fact m * Stirling (2 * n) m + P' m) / (m + 1) = -(real_of_int (int (fact m * Stirling (2 * n) m)) / real_of_int 4)" using ‹n > 0› by (auto simp: P'_def) also have "… ∈ ℤ" by (intro Ints_minus of_int_divide_in_Ints dvd) finally show ?case . next assume composite: "¬prime (m + 1)" and "m > 3" obtain a b where ab: "a * b = m + 1" "a > 1" "b > 1" using ‹m > 3› composite composite_imp_factors_nat[of "m + 1"] by auto have "a = b ⟶ a > 2" proof assume "a = b" hence "a ^ 2 > 2 ^ 2" using ‹m > 3› and ab by (auto simp: power2_eq_square) thus "a > 2" using power_less_imp_less_base by blast qed hence dvd: "(m + 1) dvd fact m" using product_dvd_fact[of a b] ab by auto have "real_of_int ((- 1) ^ m * fact m * Stirling (2 * n) m + P' m) / real (m + 1) = real_of_int ((- 1) ^ m * Stirling (2 * n) m) * (real (fact m) / (m + 1))" using composite by (auto simp: P'_def) also have "… ∈ ℤ" by (intro Ints_mult Ints_real_of_nat_divide dvd) auto finally show ?case . next assume prime: "prime (m + 1)" have "real_of_int ((-1) ^ m * fact m * int (Stirling (2 * n) m)) = (∑j≤m. (-1) ^ m * (-1) ^ (m - j) * (m choose j) * real_of_int j ^ (2 * n))" by (simp add: Stirling_closed_form sum_divide_distrib sum_distrib_left mult_ac) also have "… = real_of_int (∑j≤m. (-1) ^ j * (m choose j) * j ^ (2 * n))" unfolding of_int_sum by (intro sum.cong) (auto simp: uminus_power_if) finally have "(-1) ^ m * fact m * int (Stirling (2 * n) m) = (∑j≤m. (-1) ^ j * (m choose j) * j ^ (2 * n))" by linarith also have "… = (∑j<m+1. (-1) ^ j * (m choose j) * j ^ (2 * n))" by (intro sum.cong) auto also have "[… = (if m dvd 2 * n then - 1 else 0)] (mod (m + 1))" using vonStaudt_Clausen_lemma[of n "m + 1"] prime ‹n > 0› by simp also have "(if m dvd 2 * n then - 1 else 0) = - P' m" using prime by (simp add: P'_def) finally have "int (m + 1) dvd ((- 1) ^ m * fact m * int (Stirling (2 * n) m) + P' m)" by (simp add: cong_iff_dvd_diff) hence "real_of_int ((-1)^m * fact m * int (Stirling (2*n) m) + P' m) / of_int (int (m+1)) ∈ ℤ" by (intro of_int_divide_in_Ints) thus ?case by simp qed qed finally show ?thesis . qed subsection ‹Denominators of Bernoulli numbers› text ‹ A consequence of the Von Staudt--Clausen theorem is that the denominator of $B_{2n}$ for $n > 0$ is precisely the product of all prime numbers ‹p› such that ‹p - 1› divides $2n$. Since the denominator is obvious in all other cases, this fully characterises the denominator of Bernoulli numbers. › definition bernoulli_denom :: "nat ⇒ nat" where "bernoulli_denom n = (if n = 1 then 2 else if n = 0 ∨ odd n then 1 else ∏{p. prime p ∧ (p - 1) dvd n})" definition bernoulli_num :: "nat ⇒ int" where "bernoulli_num n = ⌊bernoulli n * bernoulli_denom n⌋" lemma finite_bernoulli_denom_set: "n > (0 :: nat) ⟹ finite {p. prime p ∧ (p - 1) dvd n}" by (rule finite_subset[of _ "{..2*n+1}"]) (auto dest!: dvd_imp_le) lemma bernoulli_denom_0 [simp]: "bernoulli_denom 0 = 1" and bernoulli_denom_1 [simp]: "bernoulli_denom 1 = 2" and bernoulli_denom_Suc_0 [simp]: "bernoulli_denom (Suc 0) = 2" and bernoulli_denom_odd [simp]: "n ≠ 1 ⟹ odd n ⟹ bernoulli_denom n = 1" and bernoulli_denom_even: "n > 0 ⟹ even n ⟹ bernoulli_denom n = ∏{p. prime p ∧ (p - 1) dvd n}" by (auto simp: bernoulli_denom_def) lemma bernoulli_denom_pos: "bernoulli_denom n > 0" by (auto simp: bernoulli_denom_def intro!: prod_pos) lemma bernoulli_denom_nonzero [simp]: "bernoulli_denom n ≠ 0" using bernoulli_denom_pos[of n] by simp lemma bernoulli_denom_code [code]: "bernoulli_denom n = (if n = 1 then 2 else if n = 0 ∨ odd n then 1 else prod_list (filter (λp. (p - 1) dvd n) (primes_upto (n + 1))))" (is "_ = ?rhs") proof (cases "even n ∧ n > 0") case True hence "?rhs = prod_list (filter (λp. (p - 1) dvd n) (primes_upto (n + 1)))" by auto also have "… = ∏(set (filter (λp. (p - 1) dvd n) (primes_upto (n + 1))))" by (subst prod.distinct_set_conv_list) auto also have "(set (filter (λp. (p - 1) dvd n) (primes_upto (n + 1)))) = {p∈{..n+1}. prime p ∧ (p - 1) dvd n}" by (auto simp: set_primes_upto) also have "… = {p. prime p ∧ (p - 1) dvd n}" using True by (auto dest: dvd_imp_le) also have "∏… = bernoulli_denom n" using True by (simp add: bernoulli_denom_even) finally show ?thesis .. qed auto corollary%important bernoulli_denom_correct: obtains a :: int where "coprime a (bernoulli_denom m)" "bernoulli m = of_int a / of_nat (bernoulli_denom m)" proof - consider "m = 0" | "m = 1" | "odd m" "m ≠ 1" | "even m" "m > 0" by auto thus ?thesis proof cases assume "m = 0" thus ?thesis by (intro that[of 1]) (auto simp: bernoulli_denom_def) next assume "m = 1" thus ?thesis by (intro that[of "-1"]) (auto simp: bernoulli_denom_def) next assume "odd m" "m ≠ 1" thus ?thesis by (intro that[of 0]) (auto simp: bernoulli_denom_def bernoulli_odd_eq_0) next assume "even m" "m > 0" define n where "n = m div 2" have [simp]: "m = 2 * n" and n: "n > 0" using ‹even m› ‹m > 0› by (auto simp: n_def intro!: Nat.gr0I) obtain a b where ab: "bernoulli (2 * n) = a / b" "coprime a (int b)" "b > 0" using Rats_int_div_natE[OF bernoulli_in_Rats] by metis define P where "P = {p. prime p ∧ (p - 1) dvd (2 * n)}" have "finite P" unfolding P_def using n by (intro finite_bernoulli_denom_set) auto from vonStaudt_Clausen[of n] obtain k where k: "bernoulli (2 * n) + (∑p∈P. 1/p) = of_int k" using ‹n > 0› by (auto simp: P_def Ints_def) define c where "c = (∑p∈P. ∏(P-{p}))" from ‹finite P› have "(∑p∈P. 1 / p) = c / ∏P" by (subst sum_inverses_conv_fraction) (auto simp: P_def prime_gt_0_nat c_def) moreover have P_nz: "prod real P > 0" using prime_gt_0_nat by (auto simp: P_def intro!: prod_pos) ultimately have eq: "bernoulli (2 * n) = (k * ∏P - c) / ∏P" using ab P_nz by (simp add: field_simps k [symmetric]) have "gcd (k * ∏P - int c) (∏P) = gcd (int c) (∏P)" by (simp add: gcd_diff_dvd_left1) also have "… = int (gcd c (∏P))" by (simp flip: gcd_int_int_eq) also have "coprime c (∏P)" unfolding c_def using ‹finite P› by (intro sum_prime_inverses_fraction_coprime) (auto simp: P_def) hence "gcd c (∏P) = 1" by simp finally have coprime: "coprime (k * ∏P - int c) (∏P)" by (simp only: coprime_iff_gcd_eq_1) have eq': "∏P = bernoulli_denom (2 * n)" using n by (simp add: bernoulli_denom_def P_def) show ?thesis by (rule that[of "k * ∏P - int c"]) (use eq eq' coprime in simp_all) qed qed lemma bernoulli_conv_num_denom: "bernoulli n = bernoulli_num n / bernoulli_denom n" (is ?th1) and coprime_bernoulli_num_denom: "coprime (bernoulli_num n) (bernoulli_denom n)" (is ?th2) proof - obtain a :: int where a: "coprime a (bernoulli_denom n)" "bernoulli n = a / bernoulli_denom n" using bernoulli_denom_correct[of n] by blast thus ?th1 by (simp add: bernoulli_num_def) with a show ?th2 by auto qed text ‹ Two obvious consequences from this are that the denominators of all odd Bernoulli numbers except for the first one are squarefree and multiples of 6: › lemma six_divides_bernoulli_denom: assumes "even n" "n > 0" shows "6 dvd bernoulli_denom n" proof - from assms have "∏{2, 3} dvd ∏{p. prime p ∧ (p - 1) dvd n}" by (intro prod_dvd_prod_subset finite_bernoulli_denom_set) auto with assms show ?thesis by (simp add: bernoulli_denom_even) qed lemma squarefree_bernoulli_denom: "squarefree (bernoulli_denom n)" by (auto intro!: squarefree_prod_coprime primes_coprime simp: bernoulli_denom_def squarefree_prime) text ‹ Furthermore, the denominator of $B_n$ divides $2(2^n - 1)$. This also gives us an upper bound on the denominators. › lemma bernoulli_denom_dvd: "bernoulli_denom n dvd (2 * (2 ^ n - 1))" proof (cases "even n ∧ n > 0") case True hence "bernoulli_denom n = ∏{p. prime p ∧ (p - 1) dvd n}" by (auto simp: bernoulli_denom_def) also have "… dvd (2 * (2 ^ n - 1))" proof (rule prime_prod_dvdI; clarify?) from True show "finite {p. prime p ∧ (p - 1) dvd n}" by (intro finite_bernoulli_denom_set) auto next fix p assume p: "prime p" "(p - 1) dvd n" show "p dvd (2 * (2 ^ n - 1))" proof (cases "p = 2") case False with p have "p > 2" using prime_gt_1_nat[of p] by force have "[2 ^ n - 1 = 1 - 1] (mod p)" using p ‹p > 2› prime_odd_nat by (intro cong_diff_nat Carmichael_divides) (auto simp: Carmichael_prime) hence "p dvd (2 ^ n - 1)" by (simp add: cong_0_iff) thus ?thesis by simp qed auto qed auto finally show ?thesis . qed (auto simp: bernoulli_denom_def) corollary bernoulli_bound: assumes "n > 0" shows "bernoulli_denom n ≤ 2 * (2 ^ n - 1)" proof - from assms have "2 ^ n > (1 :: nat)" by (intro one_less_power) auto thus ?thesis by (intro dvd_imp_le[OF bernoulli_denom_dvd]) auto qed text ‹ It can also be shown fairly easily from the von Staudt--Clausen theorem that if ‹p› is prime and ‹2p + 1› is not, then $B_{2p} \equiv \frac{1}{6}\ (\text{mod}\ 1)$ or, equivalently, the denominator of $B_{2p}$ is 6 and the numerator is of the form $6k+1$. This is the case e.\,g.\ for any primes of the form $3k+1$ or $5k+2$. › lemma bernoulli_denom_prime_nonprime: assumes "prime p" "¬prime (2 * p + 1)" shows "bernoulli (2 * p) - 1 / 6 ∈ ℤ" "[bernoulli_num (2 * p) = 1] (mod 6)" "bernoulli_denom (2 * p) = 6" proof - from assms have "p > 0" using prime_gt_0_nat by auto define P where "P = {q. prime q ∧ (q - 1) dvd (2 * p)}" have P_eq: "P = {2, 3}" proof (intro equalityI subsetI) fix q assume "q ∈ P" hence q: "prime q" "(q - 1) dvd (2 * p)" by (simp_all add: P_def) have "q - 1 ∈ {1, 2, p, 2 * p}" proof - obtain b c where bc: "b dvd 2" "c dvd p" "q - 1 = b * c" using division_decomp[OF q(2)] by auto from bc have "b ∈ {1, 2}" and "c ∈ {1, p}" using prime_nat_iff two_is_prime_nat ‹prime p› by blast+ with bc show ?thesis by auto qed hence "q ∈ {2, 3, p + 1, 2 * p + 1}" using prime_gt_0_nat[OF ‹prime q›] by force moreover have "q ≠ p + 1" proof assume [simp]: "q = p + 1" have "even q ∨ even p" by auto with ‹prime q› and ‹prime p› have "p = 2" using prime_odd_nat[of p] prime_odd_nat[of q] prime_gt_1_nat[of p] prime_gt_1_nat[of q] by force with assms show False by (simp add: cong_def) qed ultimately show "q ∈ {2, 3}" using assms ‹prime q› by auto qed (auto simp: P_def) show [simp]: "bernoulli_denom (2 * p) = 6" using ‹p > 0› P_eq by (subst bernoulli_denom_even) (auto simp: P_def) have "bernoulli (2 * p) + 5 / 6 ∈ ℤ" using ‹p > 0› P_eq vonStaudt_Clausen[of p] by (auto simp: P_def) hence "bernoulli (2 * p) + 5 / 6 - 1 ∈ ℤ" by (intro Ints_diff) auto thus "bernoulli (2 * p) - 1 / 6 ∈ ℤ" by simp then obtain a where "of_int a = bernoulli (2 * p) - 1 / 6" by (elim Ints_cases) auto hence "real_of_int a = real_of_int (bernoulli_num (2 * p) - 1) / 6" by (auto simp: bernoulli_conv_num_denom) hence "bernoulli_num (2 * p) - 1 = 6 * a" by simp thus "[bernoulli_num (2 * p) = 1] (mod 6)" by (auto simp: cong_iff_dvd_diff) qed subsection ‹Akiyama--Tanigawa algorithm› text ‹ First, we define the Akiyama--Tanigawa number triangle as shown by Kaneko~\<^cite>‹"kaneko2000"›. We define this generically, parametrised by the first row. This makes the proofs a little bit more modular. › fun gen_akiyama_tanigawa :: "(nat ⇒ real) ⇒ nat ⇒ nat ⇒ real" where "gen_akiyama_tanigawa f 0 m = f m" | "gen_akiyama_tanigawa f (Suc n) m = real (Suc m) * (gen_akiyama_tanigawa f n m - gen_akiyama_tanigawa f n (Suc m))" lemma gen_akiyama_tanigawa_0 [simp]: "gen_akiyama_tanigawa f 0 = f" by (simp add: fun_eq_iff) text ‹ The ``regular'' Akiyama--Tanigawa triangle is the one that is used for reading off Bernoulli numbers: › definition akiyama_tanigawa where "akiyama_tanigawa = gen_akiyama_tanigawa (λn. 1 / real (Suc n))" context begin private definition AT_fps :: "(nat ⇒ real) ⇒ nat ⇒ real fps" where "AT_fps f n = (fps_X - 1) * Abs_fps (gen_akiyama_tanigawa f n)" private lemma AT_fps_Suc: "AT_fps f (Suc n) = (fps_X - 1) * fps_deriv (AT_fps f n)" proof (rule fps_ext) fix m :: nat show "AT_fps f (Suc n) $ m = ((fps_X - 1) * fps_deriv (AT_fps f n)) $ m" by (cases m) (simp_all add: AT_fps_def fps_deriv_def algebra_simps) qed private lemma AT_fps_altdef: "AT_fps f n = (∑m≤n. fps_const (real (Stirling n m)) * (fps_X - 1)^m * (fps_deriv ^^ m) (AT_fps f 0))" proof - have "AT_fps f n = (fps_XD' (fps_X - 1) ^^ n) (AT_fps f 0)" by (induction n) (simp_all add: AT_fps_Suc fps_XD'_def) also have "… = (∑m≤n. fps_const (real (Stirling n m)) * (fps_X - 1) ^ m * (fps_deriv ^^ m) (AT_fps f 0))" by (rule fps_XD'_funpow_affine) simp_all finally show ?thesis . qed private lemma AT_fps_0_nth: "AT_fps f 0 $ n = (if n = 0 then -f 0 else f (n - 1) - f n)" by (simp add: AT_fps_def algebra_simps) text ‹ The following fact corresponds to Proposition 1 in Kaneko's proof: › lemma gen_akiyama_tanigawa_n_0: "gen_akiyama_tanigawa f n 0 = (∑k≤n. (- 1) ^ k * fact k * real (Stirling (Suc n) (Suc k)) * f k)" proof (cases "n = 0") case False note [simp del] = gen_akiyama_tanigawa.simps have "gen_akiyama_tanigawa f n 0 = -(AT_fps f n $ 0)" by (simp add: AT_fps_def) also have "AT_fps f n $ 0 = (∑k≤n. real (Stirling n k) * (- 1) ^ k * (fact k * AT_fps f 0 $ k))" by (subst AT_fps_altdef) (simp add: fps_sum_nth fps_nth_power_0 fps_0th_higher_deriv) also have "… = (∑k≤n. real (Stirling n k) * (- 1) ^ k * (fact k * (f (k - 1) - f k)))" using False by (intro sum.cong refl) (auto simp: Stirling_n_0 AT_fps_0_nth) also have "… = (∑k≤n. fact k * (real (Stirling n k) * (- 1) ^ k) * f (k - 1)) - (∑k≤n. fact k * (real (Stirling n k) * (- 1) ^ k) * f k)" (is "_ = sum ?f _ - ?S2</