Theory Bernoulli_FPS
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
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
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") by (simp add: sum_subtractf algebra_simps)
also from False have "sum ?f {..n} = sum ?f {0<..n}"
by (intro sum.mono_neutral_right) (auto simp: Stirling_n_0)
also have "… = sum ?f {0<..Suc n}"
by (intro sum.mono_neutral_left) auto
also have "{0<..Suc n} = {Suc 0..Suc n}" by auto
also have "sum ?f … = sum (λn. ?f (Suc n)) {0..n}"
by (subst sum.atLeast_Suc_atMost_Suc_shift) simp_all
also have "{0..n} = {..n}" by auto
also have "sum (λn. ?f (Suc n)) … - ?S2 =
(∑k≤n. -((-1)^k * fact k * real (Stirling (Suc n) (Suc k)) * f k))"
by (subst sum_subtractf [symmetric], intro sum.cong) (simp_all add: algebra_simps)
also have "-… = (∑k≤n. ((-1)^k * fact k * real (Stirling (Suc n) (Suc k)) * f k))"
by (simp add: sum_negf)
finally show ?thesis .
qed simp_all
text ‹
The following lemma states that for $A(x) := \sum_{k=0}^\infty a_{0,k} x^k$, we have
\[\sum_{n=0}^\infty a_{n,0}\frac{x^n}{n!} = e^x A(1 - e^x)\]
which correspond's to Kaneko's remark at the end of Section 2. This seems to be easier
to formalise than his actual proof of his Theorem 1, since his proof contains
an infinite sum of formal power series, and it was unclear to us how to capture this
formally.
›
lemma gen_akiyama_tanigawa_fps:
"Abs_fps (λn. gen_akiyama_tanigawa f n 0 / fact n) = fps_exp 1 * fps_compose (Abs_fps f) (1 - fps_exp 1)"
proof (rule fps_ext)
fix n :: nat
have "(fps_const (fact n) *
(fps_compose (Abs_fps (λn. gen_akiyama_tanigawa f 0 n)) (1 - fps_exp 1) * fps_exp 1)) $ n =
(∑m≤n. ∑k≤m. (1 - fps_exp 1) ^ k $ m * fact n / fact (n - m) * f k)"
unfolding fps_mult_left_const_nth
by (simp add: fps_times_def fps_compose_def gen_akiyama_tanigawa_n_0 sum_Stirling_binomial
field_simps sum_distrib_left sum_distrib_right atLeast0AtMost
del: Stirling.simps of_nat_Suc)
also have "… = (∑m≤n. ∑k≤m. (-1)^k * fact k * real (Stirling m k) * real (n choose m) * f k)"
proof (intro sum.cong refl, goal_cases)
case (1 m k)
have "(1 - fps_exp 1 :: real fps) ^ k = (-fps_exp 1 + 1 :: real fps) ^ k" by simp
also have "… = (∑i≤k. of_nat (k choose i) * (-1) ^ i * fps_exp (real i))"
by (subst binomial_ring) (simp add: atLeast0AtMost power_minus' fps_exp_power_mult mult.assoc)
also have "… = (∑i≤k. fps_const (real (k choose i) * (-1) ^ i) * fps_exp (real i))"
by (simp add: fps_const_mult [symmetric] fps_of_nat fps_const_power [symmetric]
fps_const_neg [symmetric] del: fps_const_mult fps_const_power fps_const_neg)
also have "… $ m = (∑i≤k. real (k choose i) * (- 1) ^ i * real i ^ m) / fact m"
(is "_ = ?S / _") by (simp add: fps_sum_nth sum_divide_distrib [symmetric])
also have "?S = (-1) ^ k * (∑i≤k. (-1) ^ (k - i) * real (k choose i) * real i ^ m)"
by (subst sum_distrib_left, intro sum.cong refl) (auto simp: minus_one_power_iff)
also have "(∑i≤k. (-1) ^ (k - i) * real (k choose i) * real i ^ m) =
real (Stirling m k) * fact k"
by (subst Stirling_closed_form) (simp_all add: field_simps)
finally have *: "(1 - fps_exp 1 :: real fps) ^ k $ m * fact n / fact (n - m) =
(- 1) ^ k * fact k * real (Stirling m k) * real (n choose m)"
using 1 by (simp add: binomial_fact del: of_nat_Suc)
show ?case using 1 by (subst *) simp
qed
also have "… = (∑m≤n. ∑k≤n. (- 1) ^ k * fact k *
real (Stirling m k) * real (n choose m) * f k)"
by (rule sum.cong[OF refl], rule sum.mono_neutral_left) auto
also have "… = (∑k≤n. ∑m≤n. (- 1) ^ k * fact k *
real (Stirling m k) * real (n choose m) * f k)"
by (rule sum.swap)
also have "… = gen_akiyama_tanigawa f n 0"
by (simp add: gen_akiyama_tanigawa_n_0 sum_Stirling_binomial sum_distrib_left sum_distrib_right
mult.assoc atLeast0AtMost del: Stirling.simps)
finally show "Abs_fps (λn. gen_akiyama_tanigawa f n 0 / fact n) $ n =
(fps_exp 1 * (Abs_fps f oo 1 - fps_exp 1)) $ n"
by (subst (asm) fps_mult_left_const_nth) (simp add: field_simps del: of_nat_Suc)
qed
text ‹
As Kaneko notes in his afore-mentioned remark, if we let $a_{0,k} = \frac{1}{k+1}$, we obtain
\[A(z) = \sum_{k=0}^\infty \frac{x^k}{k+1} = -\frac{\ln (1 - x)}{x}\]
and therefore
\[\sum_{n=0}^\infty a_{n,0} \frac{x^n}{n!} = \frac{x e^x}{e^x - 1} = \frac{x}{1 - e^{-x}},\]
which immediately gives us the connection to the positive Bernoulli numbers.
›
theorem bernoulli'_conv_akiyama_tanigawa: "bernoulli' n = akiyama_tanigawa n 0"
proof -
define f where "f = (λn. 1 / real (Suc n))"
note gen_akiyama_tanigawa_fps[of f]
also {
have "fps_ln 1 = fps_X * Abs_fps (λn. (-1)^n / real (Suc n))"
by (intro fps_ext) (simp del: of_nat_Suc add: fps_ln_def)
hence "fps_ln 1 / fps_X = Abs_fps (λn. (-1)^n / real (Suc n))"
by (metis fps_X_neq_zero nonzero_mult_div_cancel_left)
also have "fps_compose … (-fps_X) = Abs_fps f"
by (simp add: fps_compose_uminus' fps_eq_iff f_def)
finally have "Abs_fps f = fps_compose (fps_ln 1 / fps_X) (-fps_X)" ..
also have "fps_ln 1 / fps_X oo - fps_X oo 1 - fps_exp (1::real) = fps_ln 1 / fps_X oo fps_exp 1 - 1"
by (subst fps_compose_assoc [symmetric])
(simp_all add: fps_compose_uminus)
also have "… = (fps_ln 1 oo fps_exp 1 - 1) / (fps_exp 1 - 1)"
by (subst fps_compose_divide_distrib) auto
also have "… = fps_X / (fps_exp 1 - 1)" by (simp add: fps_ln_fps_exp_inv fps_inv_fps_exp_compose)
finally have "Abs_fps f oo 1 - fps_exp 1 = fps_X / (fps_exp 1 - 1)" .
}
also have "fps_exp (1::real) - 1 = (1 - fps_exp (-1)) * fps_exp 1"
by (simp add: algebra_simps fps_exp_add_mult [symmetric])
also have "fps_exp 1 * (fps_X / …) = bernoulli'_fps" unfolding bernoulli'_fps_def
by (subst dvd_div_mult2_eq) (auto simp: fps_dvd_iff intro!: subdegree_leI)
finally have "Abs_fps (λn. gen_akiyama_tanigawa f n 0 / fact n) = bernoulli'_fps" .
thus ?thesis by (simp add: fps_eq_iff akiyama_tanigawa_def f_def)
qed
theorem bernoulli_conv_akiyama_tanigawa:
"bernoulli n = akiyama_tanigawa n 0 - (if n = 1 then 1 else 0)"
using bernoulli'_conv_akiyama_tanigawa[of n] by (auto simp: bernoulli_conv_bernoulli')
end
end
subsection ‹Efficient code›
text ‹
We can now compute parts of the Akiyama--Tanigawa (and thereby Bernoulli numbers)
with reasonable efficiency but iterating the recurrence row by row. We essentially
start with some finite prefix of the zeroth row, say of length $n$, and then apply
the recurrence one to get a prefix of the first row of length $n - 1$ etc.
›
fun akiyama_tanigawa_step_aux :: "nat ⇒ real list ⇒ real list" where
"akiyama_tanigawa_step_aux m (x # y # xs) =
real m * (x - y) # akiyama_tanigawa_step_aux (Suc m) (y # xs)"
| "akiyama_tanigawa_step_aux m xs = []"
lemma length_akiyama_tanigawa_step_aux [simp]:
"length (akiyama_tanigawa_step_aux m xs) = length xs - 1"
by (induction m xs rule: akiyama_tanigawa_step_aux.induct) simp_all
lemma akiyama_tanigawa_step_aux_eq_Nil_iff [simp]:
"akiyama_tanigawa_step_aux m xs = [] ⟷ length xs < 2"
by (subst length_0_conv [symmetric]) auto
lemma nth_akiyama_tanigawa_step_aux:
"n < length xs - 1 ⟹
akiyama_tanigawa_step_aux m xs ! n = real (m + n) * (xs ! n - xs ! Suc n)"
proof (induction m xs arbitrary: n rule: akiyama_tanigawa_step_aux.induct)
case (1 m x y xs n)
thus ?case by (cases n) auto
qed auto
definition gen_akiyama_tanigawa_row where
"gen_akiyama_tanigawa_row f n l u = map (gen_akiyama_tanigawa f n) [l..<u]"
lemma length_gen_akiyama_tanigawa_row [simp]: "length (gen_akiyama_tanigawa_row f n l u) = u - l"
by (simp add: gen_akiyama_tanigawa_row_def)
lemma gen_akiyama_tanigawa_row_eq_Nil_iff [simp]:
"gen_akiyama_tanigawa_row f n l u = [] ⟷ l ≥ u"
by (auto simp add: gen_akiyama_tanigawa_row_def)
lemma nth_gen_akiyama_tanigawa_row:
"i < u - l ⟹ gen_akiyama_tanigawa_row f n l u ! i = gen_akiyama_tanigawa f n (i + l)"
by (simp add: gen_akiyama_tanigawa_row_def add_ac)
lemma gen_akiyama_tanigawa_row_0 [code]:
"gen_akiyama_tanigawa_row f 0 l u = map f [l..<u]"
by (simp add: gen_akiyama_tanigawa_row_def)
lemma gen_akiyama_tanigawa_row_Suc [code]:
"gen_akiyama_tanigawa_row f (Suc n) l u =
akiyama_tanigawa_step_aux (Suc l) (gen_akiyama_tanigawa_row f n l (Suc u))"
by (rule nth_equalityI) (auto simp: nth_gen_akiyama_tanigawa_row nth_akiyama_tanigawa_step_aux)
lemma gen_akiyama_tanigawa_row_numeral:
"gen_akiyama_tanigawa_row f (numeral n) l u =
akiyama_tanigawa_step_aux (Suc l) (gen_akiyama_tanigawa_row f (pred_numeral n) l (Suc u))"
by (simp only: numeral_eq_Suc gen_akiyama_tanigawa_row_Suc)
lemma gen_akiyama_tanigawa_code [code]:
"gen_akiyama_tanigawa f n k = hd (gen_akiyama_tanigawa_row f n k (Suc k))"
by (subst hd_conv_nth) (auto simp: nth_gen_akiyama_tanigawa_row length_0_conv [symmetric])
definition akiyama_tanigawa_row where
"akiyama_tanigawa_row n l u = map (akiyama_tanigawa n) [l..<u]"
lemma length_akiyama_tanigawa_row [simp]: "length (akiyama_tanigawa_row n l u) = u - l"
by (simp add: akiyama_tanigawa_row_def)
lemma akiyama_tanigawa_row_eq_Nil_iff [simp]:
"akiyama_tanigawa_row n l u = [] ⟷ l ≥ u"
by (auto simp add: akiyama_tanigawa_row_def)
lemma nth_akiyama_tanigawa_row:
"i < u - l ⟹ akiyama_tanigawa_row n l u ! i = akiyama_tanigawa n (i + l)"
by (simp add: akiyama_tanigawa_row_def add_ac)
lemma akiyama_tanigawa_row_0 [code]:
"akiyama_tanigawa_row 0 l u = map (λn. inverse (real (Suc n))) [l..<u]"
by (simp add: akiyama_tanigawa_row_def akiyama_tanigawa_def divide_simps)
lemma akiyama_tanigawa_row_Suc [code]:
"akiyama_tanigawa_row (Suc n) l u =
akiyama_tanigawa_step_aux (Suc l) (akiyama_tanigawa_row n l (Suc u))"
by (rule nth_equalityI) (auto simp: nth_akiyama_tanigawa_row
nth_akiyama_tanigawa_step_aux akiyama_tanigawa_def)
lemma akiyama_tanigawa_row_numeral:
"akiyama_tanigawa_row (numeral n) l u =
akiyama_tanigawa_step_aux (Suc l) (akiyama_tanigawa_row (pred_numeral n) l (Suc u))"
by (simp only: numeral_eq_Suc akiyama_tanigawa_row_Suc)
lemma akiyama_tanigawa_code [code]:
"akiyama_tanigawa n k = hd (akiyama_tanigawa_row n k (Suc k))"
by (subst hd_conv_nth) (auto simp: nth_akiyama_tanigawa_row length_0_conv [symmetric])
lemma bernoulli_code [code]:
"bernoulli n =
(if n = 0 then 1 else if n = 1 then -1/2 else if odd n then 0 else akiyama_tanigawa n 0)"
proof (cases "n = 0 ∨ n = 1 ∨ odd n")
case False
thus ?thesis by (auto simp add: bernoulli_conv_akiyama_tanigawa)
qed (auto simp: bernoulli_odd_eq_0)
lemma bernoulli'_code [code]:
"bernoulli' n =
(if n = 0 then 1 else if n = 1 then 1/2 else if odd n then 0 else akiyama_tanigawa n 0)"
by (simp add: bernoulli'_def bernoulli_code)
text ‹
Evaluation with the simplifier is much slower than by reflection, but can still be done
with much better efficiency than before:
›
lemmas eval_bernoulli =
akiyama_tanigawa_code akiyama_tanigawa_row_numeral
numeral_2_eq_2 [symmetric] akiyama_tanigawa_row_Suc upt_conv_Cons
akiyama_tanigawa_row_0 bernoulli_code[of "numeral n" for n]
lemmas eval_bernoulli' = eval_bernoulli bernoulli'_code[of "numeral n" for n]
lemmas eval_bernpoly =
bernpoly_def atMost_nat_numeral power_eq_if binomial_fact fact_numeral eval_bernoulli
lemma bernoulli_upto_20 [simp]:
"bernoulli 2 = 1 / 6"
"bernoulli 4 = -(1 / 30)"
"bernoulli 6 = 1 / 42"
"bernoulli 8 = - (1 / 30)"
"bernoulli 10 = 5 / 66"
"bernoulli 12 = - (691 / 2730)"
"bernoulli 14 = 7 / 6"
"bernoulli 16 = -(3617 / 510)"
"bernoulli 18 = 43867 / 798"
"bernoulli 20 = -(174611 / 330)"
by (simp_all add: eval_bernoulli)
lemma bernoulli'_upto_20 [simp]:
"bernoulli' 2 = 1 / 6"
"bernoulli' 4 = -(1 / 30)"
"bernoulli' 6 = 1 / 42"
"bernoulli' 8 = - (1 / 30)"
"bernoulli' 10 = 5 / 66"
"bernoulli' 12 = - (691 / 2730)"
"bernoulli' 14 = 7 / 6"
"bernoulli' 16 = -(3617 / 510)"
"bernoulli' 18 = 43867 / 798"
"bernoulli' 20 = -(174611 / 330)"
by (simp_all add: bernoulli'_def)
end