13 Euler-Maclaurinの和公式
黒木玄2018-07-04～2019-04-03, 2023-06-22
Copyright 2018,2019,2023 Gen Kuroki
License: MIT https://opensource.org/licenses/MIT (https://opensource.org/licenses/MIT)
Repository: https://github.com/genkuroki/Calculus (https://github.com/genkuroki/Calculus)
このファイルは次の場所できれいに閲覧できる:
http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/13%20Euler-Maclaurin%20summation%20formula.ipynb(http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/13%20Euler-Maclaurin%20summation%20formula.ipynb)
https://genkuroki.github.io/documents/Calculus/13%20Euler-Maclaurin%20summation%20formula.pdf(https://genkuroki.github.io/documents/Calculus/13%20Euler-Maclaurin%20summation%20formula.pdf)
このファイルは Julia Box (https://juliabox.com) で利用できる.
自分のパソコンにJulia言語 (https://julialang.org/)をインストールしたい場合にはWindowsへのJulia言語のインストール (http://nbviewer.jupyter.org/gist/genkuroki/81de23edcae631a995e19a2ecf946a4f)
Julia v1.1.0 の Windows 8.1 へのインストール(https://nbviewer.jupyter.org/github/genkuroki/msfd28/blob/master/install.ipynb)
を参照せよ. 前者は古く, 後者の方が新しい.
論理的に完璧な説明をするつもりはない. 細部のいい加減な部分は自分で訂正・修正せよ.

目次
 1 Bernoulli多項式
1.1 Bernoulli多項式の定義
1.2 Bernoulli多項式の基本性質
1.3 べき乗和
1.4 Bernoulli数の計算法
1.5 周期的Bernoulli多項式のFourier級数展開
 2 Euler-Maclaurinの和公式
2.1 Euler-Maclaurinの和公式の導出
2.2 Euler-Maclaurinの和公式の2つの解釈
2.3 Euler-Maclaurinの和公式の形式的導出
2.4 べき乗和(再)2.5 Stirlingの近似公式
2.6 Poissonの和公式とEuler-Maclaurinの和公式の関係
2.7 台形公式とPoissonの和公式の関係
 3 ゼータ函数への応用
3.1 解析接続
3.2 ζ(2)の近似計算
3.3 s = 1でのζ(s)の定数項がEuler定数になること
3.4 負の整数におけるゼータ函数の特殊値の計算
3.5 発散級数の有限部分と ζ(-r) の関係

In [1]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

using Base.MathConstants
using Base64
using Printf
using Statistics
const e = ℯ
endof(a) = lastindex(a)linspace(start, stop, length) = range(start, stop, length=length)​
using Plots#gr(); ENV["PLOTS_TEST"] = "true"#clibrary(:colorcet)#clibrary(:misc)
default(fmt=:png)​function pngplot(P...; kwargs...)sleep(0.1)
pngfile = tempname() * ".png"
savefig(plot(P...; kwargs...), pngfile)showimg("image/png", pngfile)end
pngplot(; kwargs...) = pngplot(plot!(; kwargs...))​
showimg(mime, fn) = open(fn) do fbase64 = base64encode(f)display("text/html", """<img src="data:$mime;base64,$base64">""")
end​
using SymPy#sympy.init_printing(order="lex") # default#sympy.init_printing(order="rev-lex")​
using SpecialFunctions
using QuadGKIn [2]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

# Override the Base.show definition of SymPy.jl:# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/​@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_fra
end@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
function toeqnarray(x::Vector{Sym})a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")"""\\left[ \\begin{array}{r}$a\\end{array} \\right]"""end
function toeqnarray(x::AbstractArray{Sym,2})sz = size(x)a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz"\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
end
print(io, as_markdown(toeqnarray(x)))
end

1 Bernoulli多項式
1.1 Bernoulli多項式の定義
定義(Bernoulli多項式): Bernoulli多項式**によって定義する.

□

𝐵𝑛 (𝑥) (𝑛 = 0,1,2,…)を
𝑧𝑒𝑧𝑥 = ∞ 𝐵𝑛 (𝑥) 𝑧𝑛
𝑒𝑧 − 1 ∑ 𝑛!

1.2 Bernoulli多項式の基本性質
一般化Bernoulli多項式の基本性質: Bernoulli多項式

𝐵0 (𝑥) = 1.
1
(2)
∫0 𝐵𝑛 (𝑥)𝑑𝑥 = 𝛿𝑛,0 .
(1)

𝑛=0

𝐵𝑛 (𝑥) は以下の性質を満たしている:

(3)

𝑛
𝑛
𝐵𝑛 (𝑥 + ℎ) = ∑ (𝑘𝑛)𝐵𝑛−𝑘(𝑥)ℎ𝑘 = ∑ (𝑘𝑛)𝐵𝑘(𝑥)ℎ𝑛−𝑘 .

𝑘=0(4) 𝐵′𝑛 (𝑥) = 𝑛𝐵𝑛−1 (𝑥) .
𝑛−1 .
(5) 𝐵𝑛 (𝑥 + 1) = 𝐵𝑛 (𝑥) + 𝑛𝑥
(6) 𝐵𝑛 (1 − 𝑥) = (−1) 𝑛 𝐵𝑛 (𝑥).
(7) 𝐵𝑛 (1) = 𝐵𝑛 (0) + 𝛿𝑛,1 となる.

𝑘=0

𝐵𝑛 (0) = 1, 𝐵𝑛 (0) = − 12 とな, 𝑛 が3以上の奇数ならば 𝐵𝑛 (0) = 0 となる.
𝑒𝑧 − 1 = 1 + 𝑂(𝑧) より, 𝑧𝑒𝑧𝑥 = 1 + 𝑂(𝑧) なので 𝐵0 (𝑥) = 1.
証明: (1) 𝑒𝑧𝑥 = 1 + 𝑂(𝑧) ,𝑧𝑒𝑧 − 1(8)(2)を示そう.

なので, これらを比較して

1

1 𝑧𝑒𝑧𝑥𝑧 1 𝑒𝑧𝑥 𝑑𝑥 = 𝑧 𝑒𝑧 − 1 = 1,𝑑𝑥=
∫0 𝑒𝑧 − 1 𝑒𝑧 − 1 ∫0𝑒𝑧 − 1 𝑧∞ 𝑧𝑛 1
1 𝑧𝑒𝑧𝑥𝐵𝑛 (𝑥)𝑑𝑥∫0 𝑒𝑧 − 1 𝑑𝑥 = ∑𝑛=0 𝑛! ∫0∫0 𝐵𝑛 (𝑥)𝑑𝑥 = 𝛿𝑛,0 .

(3) 二項定理より,

𝑛 𝑛
1
𝑛−𝑘
𝑛 𝑑𝑦 =(𝑥+𝑦)𝑥∑ (𝑘) ∫0 𝑦𝑘 𝑑𝑦.
∫0
1

𝑘=0

ゆえに,

𝑥 の函数を 𝑥 の函数に移す線形写像(前方移動平均)
1𝑓(𝑥) ↦ ∫0 𝑓(𝑥 + 𝑦)𝑑𝑦は多項式を多項式に移し, 最高次の係数が1の多項式を最高次の係数が1の同次の多項式に移す. これより, 線形写像1𝑓(𝑥) ↦ ∫0 𝑓(𝑥 + 𝑦)𝑑𝑦 は多項式どうしの一対一対応を与える線形写像になっていることがわかる. そして,
∞ ∫ 1 𝐵𝑛 (𝑥 + 𝑦)𝑑𝑦
1 𝑧𝑒𝑧(𝑥+𝑦)
0
𝑑𝑥
=
𝑧𝑛 ,
∫0 𝑒𝑧 − 1 ∑𝑛!
𝑛=0
1 𝑧𝑒𝑧(𝑥+𝑦)𝑧𝑒𝑧𝑥 1 𝑒𝑧𝑦 𝑑𝑦 = 𝑧𝑒𝑧𝑥 𝑒𝑧 − 1 = 𝑒𝑧𝑥 = ∞ 𝑥𝑛 𝑧𝑛
𝑑𝑥
=
∑ 𝑛!
∫0 𝑒𝑧 − 1𝑒𝑧 − 1 ∫0𝑒𝑧 − 1 𝑧𝑛=0なので, これらを比較して,

1∫0 𝐵𝑛 (𝑥 + 𝑦)𝑑𝑦 = 𝑥𝑛が成立することがわかる. ゆえに,

𝑛 𝑛
1 𝑛 𝑛
𝑛−𝑘 ℎ𝑘 =
𝑛=
𝑘
(𝑥
+
ℎ
+
𝑦)𝑑𝑦
=
(𝑥
+
ℎ
𝐵
)
𝑥
𝑛
∑ (𝑘)
∫0
∫0 ∑ (𝑘)𝐵𝑛−𝑘(𝑥 + 𝑦)ℎ 𝑑𝑦
1

𝑘=0

より

𝑘=0

𝑛
𝐵𝑛 (𝑥 + ℎ) = ∑ (𝑘𝑛)𝐵𝑛−𝑘(𝑥)ℎ𝑘.
𝑘=0ℎ の係数を見ることによって,𝐵′𝑛 (𝑥) = 𝑛𝐵𝑛−1 (𝑥).
(5) Bernoulli多項式の母函数の 𝑥 に 𝑥 + 1 を代入すると,
𝑧𝑒𝑧(𝑥+1) = 𝑧𝑒𝑧 𝑒𝑧𝑥 = 𝑧(1 + (𝑒𝑧 − 1))𝑒𝑧𝑥 = 𝑧𝑒𝑧𝑥 + 𝑧𝑒𝑧𝑥
𝑒𝑧 − 1 𝑒𝑧 − 1
𝑒𝑧 − 1𝑒𝑧 − 1(4) すぐ上の等式の右辺の𝑧 について展開して比較すれば(5)が得られる.
(6) Bernoulli多項式の母函数の 𝑥 に 1 − 𝑥 を代入すると,
𝑧𝑒𝑧(1−𝑥) = 𝑧𝑒𝑧 𝑒−𝑧𝑥 = 𝑧𝑒−𝑧𝑥 = −𝑧𝑒−𝑧𝑥
𝑒𝑧 − 1 𝑒𝑧 − 1 1 − 𝑒−𝑧 𝑒−𝑧 − 1とBernoulli多項式の母函数の 𝑧 に −𝑧 を代入したものになるので, 両辺を 𝑧 について展開して比較すれば(5)が得られる.
(7) 上の(2)と(4)より, 𝑛 が2以上のとき,
1
1
𝐵𝑛 (1) − 𝐵𝑛 (0) = ∫0 𝐵′𝑛 (𝑥)𝑑𝑥 = 𝑛 ∫0 𝐵𝑛−1 (𝑥)𝑑𝑥 = 𝑛𝛿𝑛−1,0 = 𝛿𝑛,1
ゆえに 𝑛 が2以上のとき 𝐵𝑛 (1) = 𝐵𝑛 (0) + 𝛿𝑛,1 .
(8) 次の函数が 𝑧 の偶函数で 𝑧 → 0 で 1 になることから, (6)が得られる:
𝑧 + 𝑧 = 𝑧 𝑒𝑧/2 + 𝑒−𝑧/2 . □𝑒𝑧 − 1 2 2 𝑒𝑧/2 − 𝑒−𝑧/2
注意: 𝐵𝑛 = 𝐵𝑛 (0) はBernoulli数と呼ばれている. (3)で (𝑥,ℎ) を (0,𝑥) で置き換えると, Bernoulli多項式がBernoulli数で表わさなので両辺をれることがわかる:𝑛𝐵𝑛 (𝑥) = ∑ (𝑘𝑛)𝐵𝑘𝑥𝑛−𝑘 .
𝑘=0

𝐵 (𝑥) が 𝑛 について帰納的に一意的に決まる □𝐵0 = 1, 𝐵1 = − 12 , 𝐵2 = 16 , 𝐵3 = 0, 𝐵4 = − 301上の定理の条件(1) (2) (4)によってBernoulli多項式例:なので

In [3]:

Out[3]:𝐵0 (𝑥) = 1, 𝐵1 (𝑥) = 𝑥 − 12 , 𝐵2 (𝑥) = 𝑥2 − 𝑥 + 16 ,𝐵3 (𝑥) = 𝑥3 − 32 𝑥2 + 12 𝑥, 𝐵4 (𝑥) = 𝑥4 − 2𝑥3 + 𝑥2 − 301 .

1 BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
2 x = symbols("x", real=true)
3 [BernoulliPolynomial(n,x) for n in 0:10]


1 

𝑥 − 12 
2 − 𝑥 + 1 𝑥
6
2 𝑥 
3𝑥3𝑥 − 2 +2


1

4
3
2
𝑥 − 2𝑥 + 𝑥 − 30 

5 − 5𝑥4 + 5𝑥3 − 𝑥 𝑥

2 4 32 6 𝑥6 − 3𝑥5 + 5𝑥2 − 𝑥2 + 421 
7 − 7𝑥6 + 7𝑥5 − 7𝑥3 + 𝑥 𝑥
26 24 62 6 14𝑥
 𝑥8 − 4𝑥7 + 3 − 7𝑥3 + 2𝑥3 − 301 
85𝑥9 − 9𝑥2 + 6𝑥7 − 21𝑥5 + 2𝑥3 − 3𝑥10 
8
2 5 
15𝑥
3𝑥
10
9
6
4
 𝑥 − 5𝑥 + 2 − 7𝑥 + 5𝑥 − 2 + 66 □

In [4]:Out[4]:
In [5]:



1
2
3
4
5

# (2) ∫_0^1 B_n(x) dx = δ_{n0}​BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
x = symbols("x", real=true)[integrate(BernoulliPolynomial(n,x), (x,0,1)) for n = 0:10]'

1
2
3
4
5
6
7

# (3) B_n(x+h) = Σ_{k=0}^n binom(n,k) B_{n-k}(x) h^k​BernoulliNumber(n) = sympy.bernoulli(n)BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
BinomCoeff(n,k) = sympy.binomial_coefficients_list(n)[k+1]x, h = symbols("x h", real=true)[BernoulliPolynomial(n,x) == sum(k->BinomCoeff(n,k)*BernoulliNumber(k)*x^(n-k), 0:n) for[1 0 0 0 0 0 0 0 0 0 0]Out[5]: 1×11 adjoint(::Vector{Bool}) with eltype Bool:
1 1 1 1 1 1 1 1 1 1 1
In [6]:



1
2
3
4
5

# (4) B_n'(x) = n B_{n-1}(x)​BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
x = symbols("x", real=true)[diff(BernoulliPolynomial(n,x), x) == n*BernoulliPolynomial(n-1,x) for n = 1:10]'

Out[6]: 1×10 adjoint(::Vector{Bool}) with eltype Bool:
1 1 1 1 1 1 1 1 1 1
In [7]:

Out[7]:

In [8]:



1
2
3
4
5

# (5) B_n(x+1) = B_n(x) + n x^{n-1}​BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
x = symbols("x", real=true)[simplify(BernoulliPolynomial(n,x+1) - BernoulliPolynomial(n,x)) for n in 0:10]

1
2
3
4
5

# (6) B_n(1-x) = (-1)^n B_n(x)​BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
x = symbols("x", real=true)[expand(BernoulliPolynomial(n,1-x)) == (-1)^n*BernoulliPolynomial(n,x) for n in 0:10]'

 0 
 1 
 3𝑥2𝑥2  4𝑥3  5𝑥4  6𝑥56  7𝑥
8𝑥7 
 9𝑥8  10𝑥9 Out[8]: 1×11 adjoint(::Vector{Bool}) with eltype Bool:
1 1 1 1 1 1 1 1 1 1 1
In [9]:

Out[9]:

1
2
3
4
5# (7) B_n(1) = B_n(0) + δ_{n1}​BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
x = symbols("x", real=true)[expand(BernoulliPolynomial(n,1)) - BernoulliPolynomial(n,0) for n in 0:10]'[0 1 0 0 0 0 0 0 0 0 0]

In [10]:



1
2
3
4

# (8) B_n = B_n(0) は n が3以上の奇数ならば0になる.
​
BernoulliNumber(n) = sympy.bernoulli(n)[(n, BernoulliNumber(n)) for n in 0:10]

Out[10]: 11-element Vector{Tuple{Int64, Sym}}:(0, 1)(1, -1/2)(2, 1/6)(3, 0)(4, -1/30)(5, 0)(6, 1/42)(7, 0)(8, -1/30)(9, 0)(10, 5/66)

1.3 べき乗和

𝑚 は正の整数であるする. Bernoulli多項式について,
𝐵𝑚+1(𝑥 + 1) − 𝐵𝑚+1(𝑥) = (𝑚 + 1)𝑥𝑚 , i.e. 𝑥𝑚 = 𝐵𝑚+1(𝑥 +𝑚1)+−1 𝐵𝑚+1(𝑥)が成立しているので, これを 𝑥 = 0,1,…,𝑛 について足し上げると,𝑛
𝐵𝑚+1(𝑛 + 1) − 𝐵𝑚+1 . □∑ 𝑗𝑚 =
𝑚+1
𝑗=1

In [11]:



1
2
3
4
5

PowerSum(m, n) = sum(j->j^m, 1:n)BernoulliNumber(n) = sympy.bernoulli(n)BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)
PowerSumFormula(m, n) = (BernoulliPolynomial(m+1,n+1)-BernoulliNumber(m+1))/(m+1)[(m, PowerSum(m,10), PowerSumFormula(m, 10)) for m in 1:10]

Out[11]: 10-element Vector{Tuple{Int64, Int64, Sym}}:(1, 55, 55)(2, 385, 385)
(3, 3025, 3025)(4, 25333, 25333)(5, 220825, 220825)(6, 1978405, 1978405)(7, 18080425, 18080425)(8, 167731333, 167731333)(9, 1574304985, 1574304985)(10, 14914341925, 14914341925)

1.4 Bernoulli数の計算法
Bernoulli数

𝐵𝑛 は

𝑧 = ∞ 𝐵𝑛 𝑧𝑛
𝑒𝑧 − 1 ∑𝑛=1 𝑛!

で定義される. しかし, この展開を直接計算することによって Bernoulli 数を求めるのは効率が悪い.

𝑧 → 0 の極限を取ることによって 𝐵0 = 1 であることはすぐにわかる.
次に, 𝑛 が 3 以上の奇数のとき 𝐵𝑛 = 0 となることを(再び)示そう.
𝑧 + 𝑧 = 𝑧 𝑒𝑧 + 1 = 𝑧 𝑒𝑧/2 + 𝑒−𝑧/2
𝑒𝑧 − 1 2 2 𝑒𝑧 − 1 2 𝑒𝑧/2 − 𝑒−𝑧/2
より, 左辺は偶函数になるので, その展開の奇数次の項は消える. このことから, 𝐵1 = −1/2 でかつ, 0 = 𝐵3 = 𝐵5 = 𝐵7 = ⋯まず, 左辺のであることもわかる.

𝑧𝑒𝑧 = ∞ 𝑧𝑗 𝐵𝑘𝑧𝑘 = ∞ 𝑛 𝑛 𝐵𝑘 𝑧𝑛∑ 𝑗! 𝑘! ∑𝑒𝑧 − 1 𝑗,𝑘=0
𝑛=0 (∑
𝑘=0 (𝑘) ) 𝑛!
でかつ𝑧𝑒𝑧 = 𝑧 + 𝑧 = ∞ (𝐵𝑛 + 𝛿𝑛1 ) 𝑧𝑛∑𝑛!
𝑒𝑧 − 1 𝑒𝑧 − 1𝑛=0なので, これらを比較すると

𝑛 𝐵𝑘 = 𝛿𝑛1 .
∑
(
𝑘
𝑘=0 )
𝑛−1

ゆえに,

𝑛 を 𝑛 + 1 で置き換え, 𝑛 ≧ 1 とし, 𝐵𝑛 を他で表わす式に書き直すと
𝑛−1
𝐵𝑛 = − 𝑛 +1 1 ∑ (𝑛 +𝑘 1)𝐵𝑘 (𝑛 ≧ 1).
𝑘=0

𝐵𝑛 を求めることができる. 𝐵0 = 1, 𝐵1 = −1/2, 0 = 𝐵3 = 𝐵5 = 𝐵7 = ⋯ であることを使うと,
𝑚−1
𝐵2𝑚 = − 2𝑚1+ 1 (1 − 2𝑚2+ 1 + ∑ (2𝑚2𝑘+ 1)𝐵2𝑘) .
問題: 上の方ではSymPyにおけるBernoulli数の函数を利用した. Bernoulli数を計算するためのプログラムを自分で書け. □解答例: 次のセルの通り. □これを使えば帰納的に

In [12]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

# binomial coefficient: binom(n,k) = n(n-1)・(n-k+1)/k!
#mydiv(a, b) = a / b
mydiv(a::Integer, b::Integer) = a ÷ bfunction binom(n, k)
k < 0 && return zero(n)k == 0 && return one(n)b = one(n)for j in 1:k
b = mydiv(b*(n-k+j), j)
end
b
end
@show binom(Rational(big"100")/3, 30)​# Bernoulli numbers: B(n) = Bernoulli[n+1] = B_n#struct Bernoulli{T}
B::Array{T,1}end
function Bernoulli(; maxn=200)
B = zeros(Rational{BigInt},maxn+1)
B[1] = 1# B_0
B[2] = -1//2 # B_1
for n in big"2":2:maxn+1
B[n+1] = -(1//(n+1))*sum(j->binom(n+1,j)*B[j+1], 0:n-1)# B_n = -(1/(n+1)) Σ_{j=0}^{n-1} binom(n+1,j)*B_j
end
Bernoulli(B)end(B::Bernoulli)(n) = B.B[n+1]​maxn = 200@time B = Bernoulli(maxn=maxn) # B_n を B_{maxn} まで計算
BB(n) = float(B(n)) # B(n) = B_n である. BB(n)はその浮動小数点版​# SymPyのBernoulli数と比較して正しく計算できているかどうかを確認#
BernoulliNumber(n) = sympy.bernoulli(n)@show B_eq_B = [B(n) == BernoulliNumber(n) for n in 0:maxn]
println()@show all(B_eq_B)​
maxnprint = 30
println()for n in [0; 1; 2:2:maxnprint]
println("B($n) = ", B(n))
endprintln()for n in [0; 1; 2:2:maxnprint]
println("BB($n) = ", BB(n))
endbinom(Rational(big"100") / 3, 30) = 11240781188817808072725280//984770902183611232881
1.173728 seconds (12.36 M allocations: 261.706 MiB, 19.32% gc time, 19.96% compilation time)
B_eq_B = [B(n) == BernoulliNumber(n) for n = 0:maxn] = Bool[1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
all(B_eq_B) = true
B(0) = 1//1
B(1) = -1//2
B(2) = 1//6
B(4) = -1//30
B(6) = 1//42
B(8) = -1//30
B(10) = 5//66
B(12) = -691//2730
B(14) = 7//6
B(16) = -3617//510
B(18) = 43867//798
B(20) = -174611//330
B(22) = 854513//138
B(24) = -236364091//2730
B(26) = 8553103//6
B(28) = -23749461029//870
B(30) = 8615841276005//14322
BB(0) = 1.0
BB(1) = -0.5
BB(2) = 0.1666666666666666666666666666666666666666666666666666666666666666666666666666674
BB(4) = -0.03333333333333333333333333333333333333333333333333333333333333333333333333333359
BB(6) = 0.02380952380952380952380952380952380952380952380952380952380952380952380952380947
BB(8) = -0.03333333333333333333333333333333333333333333333333333333333333333333333333333359
BB(10) = 0.0757575757575757575757575757575757575757575757575757575757575757575757575757578
BB(12) = -0.2531135531135531135531135531135531135531135531135531135531135531135531135531131
BB(14) = 1.166666666666666666666666666666666666666666666666666666666666666666666666666661
BB(16) = -7.092156862745098039215686274509803921568627450980392156862745098039215686274513
BB(18) = 54.97117794486215538847117794486215538847117794486215538847117794486215538847111
BB(20) = -529.1242424242424242424242424242424242424242424242424242424242424242424242424247
BB(22) = 6192.123188405797101449275362318840579710144927536231884057971014492753623188388
BB(24) = -86580.25311355311355311355311355311355311355311355311355311355311355311355311313
BB(26) = 1.425517166666666666666666666666666666666666666666666666666666666666666666666661e+
06
BB(28) = -2.729823106781609195402298850574712643678160919540229885057471264367816091954032e+07
BB(30) = 6.015808739006423683843038681748359167714006423683843038681748359167714006423638e+
08

1.5 周期的Bernoulli多項式のFourier級数展開

𝐵˜𝑘(𝑥) = 𝐵𝑘(𝑥 − ⌊𝑥⌋) を周期的Bernoulli多項式と呼ぶことにする. 周期的Bernoulli多項式は 𝐵˜𝑘(𝑥 + 1) = 𝐵˜𝑘(𝑥) を満たしている.

𝑧𝑒𝑧(𝑥−⌊𝑥⌋) の 𝑥 の函数としての Fourier係数 𝑎𝑛 (𝑧) は次のように求まる:
𝑒𝑧 − 1𝑒𝑧 − 1 𝑎𝑛 (𝑧) = 1 𝑒𝑧𝑥 𝑒−2𝜋𝑖𝑛𝑥 𝑑𝑥 = 𝑒(𝑧−2𝜋𝑖𝑛)𝑥 𝑥=1 = 𝑒𝑧 − 1 , 𝑎𝑛 (𝑧) = 𝑧 .
∫0[ 𝑧 − 2𝜋𝑖𝑛 ]𝑥=0 𝑧 − 2𝜋𝑖𝑛𝑧 − 2𝜋𝑖𝑛𝑧
ゆえに 𝑎0 (𝑧) = 1 であり, 𝑛 ≠ 0 のとき∞ 𝑧𝑘
𝑎𝑛 (𝑧) = − ∑ (2𝜋𝑖𝑛)𝑘𝑘=1
これより, ˜𝐵 𝑘(𝑥) のFourier係数 𝑎𝑘,𝑛 は, 𝑎0,𝑛 = 𝛿𝑛,0 , 𝑎𝑘,0 = 𝛿𝑘,0 を満たし, 𝑘 ≠ 0, 𝑛 ≧ 1 のとき
𝑘!
𝑎𝑘,𝑛 = − (2𝜋𝑖𝑛)𝑘
周期的Bernoulli多項式の母函数

となることがわかる. したがって, Fourier級数論より,
実数 について次が成立することがわかる:

𝑥

すなわち,

𝑒2𝜋𝑖𝑛𝑥 .
𝐵˜𝑘(𝑥) = 𝐵𝑘(𝑥 − ⌊𝑥⌋) = −𝑘! ∑ (2𝜋𝑖𝑛)𝑘𝑘 = 1,2,3,… について∞𝐵˜2𝑘−1(𝑥) = (−1)𝑘2(2𝑘 − 1)! ∑ sin(2𝜋𝑛𝑥)(2𝜋𝑛)2𝑘−1 ,𝑛=1

このことから,In [13]:

Out[13]:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

𝑘 = 1 のときは整数ではない実数 𝑥 について, 𝑘 ≧ 2 の場合にはすべての
𝑛≠0∞
𝐵˜2𝑘(𝑥) = (−1)𝑘−12(2𝑘)! ∑ cos(2𝜋𝑛𝑥)(2𝜋𝑛)2𝑘 .
𝑛=1

𝑘 が大きいとき(実際には 𝑘 = 5,6 程度ですでに), 周期的Bernoulli多項式は 𝑛 = 1 の項だけで𝑘−12(2𝑘)! cos(2𝜋𝑥)
˜
,
(𝑥)
≈
(−1
𝐵˜2𝑘−1(𝑥) ≈ (−1)𝑘2(2𝑘 − 1)! sin(2𝜋𝑥)
𝐵
)
2𝑘
2𝑘−1
(2𝜋)
(2𝜋)2𝑘BBB = Bernoulli(Float64.(B.B)) # Float64 Bernoulli numbers
BP(k,x) = sum(j->binom(k,j)*BBB(k-j)*x^j, 0:k) # Float64 Bernoulli polynomial
PBP(k,x) = BP(k, x - floor(x)) # periodic Bernoulli polynomial​# partial sum of Fourier series of periodic Bernoulli polynomial
function PSFS(k, N, x)
k == 0 && return zero(x)if isodd(k)return (-1)^((k+1)÷2)*2*factorial(k)*sum(n->sin(2π*n*x)/(2π*n)^k, 1:N)
else
return (-1)^(k÷2-1)*2*factorial(k)*sum(n->cos(2π*n*x)/(2π*n)^k, 1:N)
end
end
​
PP = []
x = -1.0:0.001:0.999
for (k,N) in [(1,20), (2,10), (3,3), (4,2), (5,1), (6,1)]
y = PBP.(k,x)
z = PSFS.(k, N, x)ymin = 1.2*minimum(y)ymax = 2.7*maximum(y)P = plot(legend=:topleft, size=(400, 250), ylim=(ymin, ymax))
plot!(x, y, label="B_$k(x-[x])")plot!(x, z, label="partial sum of Fourier series (N=$N)")
push!(PP, P)
end​plot(PP[1:2]..., size=(750, 280))

In [14]:

1 plot(PP[3:4]..., size=(750, 280))

Out[14]:

In [15]:

1 plot(PP[5:6]..., size=(750, 280))

Out[15]:

2 Euler-Maclaurinの和公式
2.1 Euler-Maclaurinの和公式の導出
Bernoulli多項式

𝐵𝑛 (𝑥) とBernoulli数 𝐵𝑛 について
𝐵0 (𝑥) = 1, 𝑑𝑥𝑑 𝐵𝑛𝑛!(𝑥) = 𝐵(𝑛𝑛−1− (𝑥)1)! ,𝐵1 (0) = − 12 , 𝐵1 (1) = 12 ,𝐵𝑛 (1) = 𝐵𝑛 (0) = 𝐵𝑛 (𝑛 = 0,2,3,4,5,…)
𝐵2𝑗+1 = 0 (𝑗 = 1,2,3,…)が成立している. 以下ではしばらくのあいだこれらの条件しか使わない.
部分積分を繰り返すことによって,

1

1∫0 𝑓(𝑥)𝑑𝑥 = ∫0 𝐵0 (𝑥)𝑓(𝑥)𝑑𝑥1= [𝐵1 (𝑥)𝑓(𝑥)]10 − ∫0 𝐵1 (𝑥)𝑓 ′ (𝑥)𝑑𝑥1= [𝐵1 (𝑥)𝑓(𝑥)]10 − 12 [𝐵2 (𝑥)𝑓 ′ (𝑥)]10 + ∫0 𝐵22(𝑥) 𝑓 ″ (𝑥)𝑑𝑥1= [𝐵1 (𝑥)𝑓(𝑥)]10 − 12 [𝐵2 (𝑥)𝑓 ′ (𝑥)]10 + 3!1 [𝐵3 (𝑥)𝑓 ″ (𝑥)]10 − ∫0 𝐵33!(𝑥) 𝑓 ‴ (𝑥)𝑑𝑥= ⋯⋯⋯⋯⋯𝑛1𝑘−1= ∑ (−1)𝑘! [𝐵𝑘(𝑥)𝑓 (𝑘−1) (𝑥)]10 + (−1)𝑛 ∫0 𝐵𝑛𝑛!(𝑥) 𝑓 (𝑛) (𝑥)𝑑𝑥
𝑘=1
𝑛1= 𝑓(0) +2 𝑓(1) + ∑(−1)𝑘−1 𝐵𝑘!𝑘 (𝑓 (𝑘−1) (1) − 𝑓 (𝑘−1) (0)) + (−1)𝑛 ∫0 𝐵𝑛𝑛!(𝑥) 𝑓 (𝑛) (𝑥)𝑑𝑥.
𝑘=2

𝑥

𝑥
⌊𝑥⌋ と書く. このとき, 𝑥 − ⌊𝑥⌋ は 𝑥 の「小数部分」になる. このように記号を準備し
𝑗1
𝑗+1
𝑓(𝑥)𝑑𝑥=∫𝑗∫0 𝑓(𝑥 + 𝑗)𝑑𝑥
+ 1) + 𝑛 (−1)𝑘−1 𝐵𝑘 (𝑓 (𝑘−1) (𝑗 + 1) − 𝑓 (𝑘−1) (𝑗)) + (−1)𝑛 1 𝐵𝑛 (𝑥) 𝑓 (𝑛) (𝑥 + 𝑗)𝑑𝑥
= 𝑓(𝑗) + 𝑓(𝑗
∑
∫0 𝑛!
𝑘!
2
𝑘=2
+ 1) + 𝑛 (−1)𝑘−1 𝐵𝑘 (𝑓 (𝑘−1) (𝑗 + 1) − 𝑓 (𝑘−1) (𝑗)) + (−1)𝑛 𝑗+1 𝐵𝑛 (𝑥 − ⌊𝑥⌋) 𝑓 (𝑛) (𝑥)𝑑𝑥.
= 𝑓(𝑗) + 𝑓(𝑗
∑
∫𝑗
2
𝑘!
𝑛!

実数 に対して, 以下の最大の整数をておくと, 整数 に対して,

𝑘=2

𝑎 < 𝑏 を満たす整数 𝑎,𝑏 に対して, 上の式を 𝑗 = 𝑎 から 𝑗 = 𝑏 − 1 まで足し上げると,
𝑏−1 𝑓(𝑗) + 𝑓(𝑗 + 1) 𝑛𝑏𝑏 𝐵𝑛 (𝑥 − ⌊𝑥⌋)𝑘−1 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) + (−1)𝑛𝑓(𝑥)𝑑𝑥=+(−1
𝑓 (𝑛) (𝑥)𝑑𝑥
)
∫𝑎
∑
∑
∫
𝑘!
𝑛!
2
𝑎
𝑗=𝑎
𝑘=2
𝑛
𝑏−1
= 𝑓(𝑎) +2 𝑓(𝑏) + ∑ 𝑓(𝑗) + ∑(−1)𝑘−1 𝐵𝑘!𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) − 𝑅𝑛 .
𝑗=𝑎+1
𝑘=2
ここで

𝑏
𝑅𝑛 = (−1)𝑛−1 ∫𝑎 𝐵𝑛 (𝑥 𝑛!− ⌊𝑥⌋) 𝑓 (𝑛) (𝑥)𝑑𝑥
とおいた. さらい, 𝑛 が3以上の奇数のとき 𝐵𝑛 = 0 となることを使うと,
𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑏−1 𝑓(𝑗) −𝐵
𝑓(𝑥)𝑑𝑥=∑∑ (2𝑖)!2𝑖 (𝑓 (2𝑖−1) (𝑏) − 𝑓 (2𝑖−1) (𝑎)) − 𝑅𝑛 .
∫𝑎
2
𝑗=𝑎+1

この公式をEuler-Maclaurinの和公式と呼ぶ. これはる:1≦𝑖≦𝑛/2

𝑛 が3以上の奇数のとき 𝐵𝑛 = 0 となることを使うと次のように書き直され𝑓(𝑎) + 𝑓(𝑏) + 𝑏−1 𝑓(𝑗) − 𝑛 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) − 𝑅𝑛 .
𝑓(𝑥)𝑑𝑥
=
∑
∑ 𝑘!
∫𝑎
2
𝑏

𝑗=𝑎+1

𝑘=2

𝑛 が2以上の整数の場合の結果である. 𝑛 = 1 の場合の結果は次のようにして得られる.
1
1
1
1′ 𝑓(𝑥)𝑑𝑥 = [𝑥𝑓(𝑥)]1 − 𝑥𝑓 ′ (𝑥)𝑑𝑥 = 𝑓(1) − 𝑥𝑓 ′ (𝑥)𝑑𝑥
𝑓(𝑥)𝑑𝑥
=
𝑥
0
∫0
∫0
∫0
∫0
より, 𝑎 < 𝑏 を満たす整数 𝑎,𝑏 に対して,
𝑏−1
𝑏−1 1
𝑏
𝑓(𝑥)𝑑𝑥
=
𝑓(𝑗
+
1)
−
∑
∑ ∫0 𝑥𝑓 ′ (𝑥 + 𝑗)𝑑𝑥∫𝑎
注意: 上の公式は

𝑗=𝑎
𝑏

𝑗=𝑎

𝑏

= ∑ 𝑓(𝑗) − 𝑓(𝑎) − ∫𝑎 (𝑥 − ⌊𝑥⌋)𝑓 ′ (𝑥)𝑑𝑥
𝑗=𝑎

すなわち,

𝑏

𝑏

𝑏

𝑏

𝑏−1

′
∑ 𝑓(𝑗) = ∫𝑎 𝑓(𝑥)𝑑𝑥 + 𝑓(𝑎) + ∫𝑎 (𝑥 − ⌊𝑥⌋)𝑓 (𝑥)𝑑𝑥
𝑗=𝑎1= ∫𝑎 𝑓(𝑥)𝑑𝑥 + 𝑓(𝑎) + ∑ ∫0 𝑥𝑓 ′ (𝑥 + 𝑗)𝑑𝑥

これらのシンプルな公式も結構よく使われる.

□

𝑗=𝑎

注意: Euler-Maclaurinの和公式の証明は部分積分の繰り返しであった. それはTaylorの定理の証明とほぼ同じような議論である.
Taylorの定理も以下のような部分積分の繰り返しで証明される:

𝑝𝑛 (𝑡) = (𝑥 −𝑛! 𝑡)(𝑥 − 𝑎)𝑛 なのでとおくと, 𝑝0 (𝑡) = 1 , −𝑝′𝑛 (𝑡) = 𝑝𝑛−1 (𝑡), [𝑔(𝑡)(−𝑝𝑛 (𝑡))] 𝑥𝑎 = 𝑔(𝑎)
𝑛!
𝑛

𝑥

𝑥

𝑓(𝑥) = 𝑓(𝑎) + ∫𝑎 𝑓(𝑡)𝑑𝑡 = 𝑓(𝑎) + ∫𝑎 𝑓(𝑡)𝑝0 (𝑡)𝑑𝑡𝑥= 𝑓(𝑎) + 𝑓 ′ (𝑎)(𝑥 − 𝑎) + ∫𝑎 𝑓 ′ (𝑡)𝑝1 (𝑡)𝑑𝑡𝑥= 𝑓(𝑎) + 𝑓 ′ (𝑎)(𝑥 − 𝑎) + 𝑓 ″ (𝑎) (𝑥 −2 𝑎)2 + ∫𝑎 𝑓 ″ (𝑡)𝑝2 (𝑡)𝑑𝑡𝑥= 𝑓(𝑎) + 𝑓 ′ (𝑎)(𝑥 − 𝑎) + 𝑓 ″ (𝑎) (𝑥 −2 𝑎)2 + 𝑓 ‴ (𝑎) (𝑥 −3!𝑎)3 + ∫𝑎 𝑓 ‴ (𝑡)𝑝3 (𝑡)𝑑𝑡= ⋯⋯⋯⋯⋯𝑛= ∑ 𝑓 (𝑘) (𝑎) (𝑥 −𝑘!𝑎)𝑘 + 𝑅𝑛 .
𝑘=0
ここで

𝑅𝑛 = ∫𝑎 𝑓 (𝑛) 𝑝𝑛 (𝑡)𝑑𝑡 = ∫𝑎 𝑓 (𝑛) (𝑥 −𝑛! 𝑡) 𝑑𝑡
とおいた この意味でEuler Maclaurinの和公式はTaylorの公式の「仲間」だと言える □𝑥𝑥𝑛

2.2 Euler-Maclaurinの和公式の2つの解釈
Euler-Maclaurinの和公式の解釈1:

𝑏−1 𝑓(𝑗) + 𝑓(𝑗 + 1)

は積分

𝑏

∑
𝑗=𝑎

2

= 𝑓(𝑎) +2 𝑓(𝑏) + ∑ 𝑓(𝑗)
𝑏−1

𝑗=𝑎+1

∫𝑎 𝑓(𝑥)𝑑𝑥 の近似計算に使われる台形公式である. ゆえに, Euler-Maclaurinの和公式
𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑏−1 𝑓(𝑗) −𝐵2𝑖 (𝑓 (2𝑖−1) (𝑏) − 𝑓 (2𝑖−1) (𝑎)) − 𝑅𝑛 ,𝑓(𝑥)𝑑𝑥
=
∑
∑
∫𝑎
(2𝑖)!
2
𝑗=𝑎+1
1≦𝑖≦𝑛/2𝑏 𝐵𝑛 (𝑥 − ⌊𝑥⌋)(𝑛)
𝑅𝑛 = (−1)𝑛−1 ∫𝑎𝑛! 𝑓 (𝑥)𝑑𝑥

は台形公式による積分の近似の誤差が,𝐵2𝑖 (𝑓 (2𝑖−1) (𝑏) − 𝑓 (2𝑖−1) (𝑎)) − 𝑅𝑛
− ∑ (2𝑖)!
1≦𝑖≦𝑛/2

𝑛 = 1 の場合には, 𝐵1 (𝑥) = 𝑥 − 12 なので,𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑏−1 𝑓(𝑗) − 𝑏 𝑥 − ⌊𝑥⌋ − 1 𝑓 ′ (𝑥)𝑑𝑥.
𝑓(𝑥)𝑑𝑥
=
∑
∫𝑎
∫𝑎 (
2
2)
𝑗=𝑎+1
𝑛 = 2 の場合には 𝐵2 (𝑥) = 𝑥2 − 𝑥 + 16 , 𝐵2 = 16 であり,𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑏−1 𝑓(𝑗) − 𝑓 ′ (𝑏) − 𝑓 ′ (𝑎) + 𝑏 𝐵2 (𝑥 − ⌊𝑥⌋) 𝑓 ″ (𝑥)𝑑𝑥.
𝑓(𝑥)𝑑𝑥
=
∑
∫𝑎
∫𝑎
2
12
2
𝑗=𝑎+1
となる. □になっていることを意味している. 例えば,

In [16]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# すぐ上の公式の検証​
BernoulliNumber(n) = sympy.bernoulli(n)BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)​
function EulerMaclaurinIntegral(f, a, b, n)
x = symbols("x", real=true)((f(a)+f(b))/Sym(2)+ sum(j->f(j), a+1:b-1)- sum(k -> (BernoulliNumber(k)/factorial(Sym(k))* (diff(f(x), x, k-1)(x=>b) - diff(f(x), x, k-1)(x=>a))), 2:n))end​
function EulerMaclaurinRemainder(f, a, b, n)
x = symbols("x", real=true)g = diff(f(x), x, n)(-1)^(n-1) * sum(k -> (integrate(BernoulliPolynomial(n,x)*g(x=>x+k), (x,0,1))), a:b-1)/factorial(Sym(n))
end​x = symbols("x", real=true)​[integrate(x^m, (x, 0, 10)) for m in 7:15] |> display​[
EulerMaclaurinIntegral(x->x^m, 0, 10, 5) - EulerMaclaurinRemainder(x->x^m, 0, 10, 5)
for m in 7:15
] |> display

12500000  1000000000

9
1000000000  100000000000 
11 
 250000000000

3
 10000000000000

13
 50000000000000 
7
 200000000000000

3


625000000000000 
12500000  1000000000 
9
1000000000  100000000000 
11 
 250000000000

3


 10000000000000
13
 50000000000000 
7
 200000000000000

3


625000000000000 

Euler-Maclaurinの和公式の解釈2: Euler-Maclaurinの和公式は次のように書き直される:𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) +𝐵2𝑖 (𝑓 (2𝑖−1) (𝑏) − 𝑓 (2𝑖−1) (𝑎)) + 𝑅𝑛 ,𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∑
∫
2
(2𝑖)!
𝑎
𝑗=𝑎
1≦𝑖≦𝑛/2𝑏 𝐵𝑛 (𝑥 − ⌊𝑥⌋)(𝑛)
𝑅𝑛 = (−1)𝑛−1 ∫𝑎𝑛! 𝑓 (𝑥)𝑑𝑥
これは 𝑛 が3以上の奇数のとき 𝐵𝑛 = 0 となることを使うと次のように書き直される:𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑛 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) + 𝑅𝑛 ,𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∫𝑎
∑
2
𝑗=𝑎
𝑘=2 𝑘!
𝑏
𝑅𝑛 = (−1)𝑛−1 ∫𝑎 𝐵𝑛 (𝑥 𝑛!− ⌊𝑥⌋) 𝑓 (𝑛) (𝑥)𝑑𝑥𝑏𝑏𝑓(𝑗)この等式は函数 𝑓 の整数における値の和を積分∑∫𝑎 𝑓(𝑥)𝑑𝑥 で近似したときの誤差が
𝑗=𝑎𝑓(𝑎) + 𝑓(𝑏) +𝐵2𝑖 (𝑓 (2𝑖−1) (𝑏) − 𝑓 (2𝑖−1) (𝑎)) + 𝑅𝑛
∑
(2𝑖)!
2
1≦𝑖≦𝑛/2
1 なので,になっていることを意味している. 例えば, 𝑛 = 1 の場合には, 𝐵1 (𝑥) = 𝑥 −2𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑏 𝑥 − ⌊𝑥⌋ − 1 𝑓 ′ (𝑥)𝑑𝑥.
𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∫𝑎
∫𝑎 (
2
2)
𝑗=𝑎
𝑛 = 2 の場合には 𝐵2 (𝑥) = 𝑥2 − 𝑥 + 16 , 𝐵2 = 16 であり,𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑓 ′ (𝑏) − 𝑓 ′ (𝑎) − 𝑏 𝐵2 (𝑥 − ⌊𝑥⌋) 𝑓 ″ (𝑥)𝑑𝑥.
𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∫
∫
2
12
2

In [17]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

# すぐ上の公式を検証​
PowerSum(m, n) = sum(j->j^m, 1:n)BernoulliNumber(n) = sympy.bernoulli(n)BernoulliPolynomial(n,x) = sympy.bernoulli(n,x)​
function EulerMaclaurinSum(f, a, b, n)
x = symbols("x", real=true)(
integrate(f(x), (x, a, b))+ (f(a)+f(b))/Sym(2)+ sum(k -> (BernoulliNumber(k)/factorial(Sym(k))* (diff(f(x), x, k-1)(x=>b) - diff(f(x), x, k-1)(x=>a))), 2:n))end​
function EulerMaclaurinRemainder(f, a, b, n)
x = symbols("x", real=true)g = diff(f(x), x, n)(-1)^(n-1) * sum(k -> (integrate(BernoulliPolynomial(n,x)*g(x=>x+k), (x,0,1))), a:b-1)/factorial(Sym(n))
end​[PowerSum(m, 10) for m in 1:10] |> display​[EulerMaclaurinSum(x->x^m, 1, 10, m+1) for m in 1:10] |> display​[
EulerMaclaurinSum(x->x^m, 1, 10, m-1) + EulerMaclaurinRemainder(x->x^m, 1, 10, m-1)
for m in 3:10
] |> display​[
EulerMaclaurinSum(x->x^m, 1, 10, m-2) + EulerMaclaurinRemainder(x->x^m, 1, 10, m-2)
for m in 4:10
] |> display​[
EulerMaclaurinSum(x->x^m, 1, 10, m-3) + EulerMaclaurinRemainder(x->x^m, 1, 10, m-3)
for m in 5:10
] |> display

10-element Vector{Int64}:
55
385
3025
25333
220825
1978405
18080425
167731333
1574304985
14914341925


55 
385 

3025 

25333
 220825  1978405 
18080425  167731333
 1574304985  14914341925 3025 
25333 
220825
 1978405  18080425  167731333  1574304985 
14914341925
25333  220825  1978405 
18080425  167731333
 1574304985  14914341925  220825  1978405 
18080425  167731333
 1574304985  14914341925 
2.3 Euler-Maclaurinの和公式の形式的導出
函数

𝑓(𝑥) に対して, ある函数 𝐹(𝑥) で𝐹(𝑥 + 1) − 𝐹(𝑥) = 𝑓(𝑥 + ℎ)
∂ とおくと, 形式的にその条件はという条件を満たすものを求める問題を考える. そのとき, 𝐷 =∂𝑥(𝑒𝐷 − 1)𝐹(𝑥) = 𝑒ℎ𝐷 𝑓(𝑥) = 𝐷𝑒ℎ𝐷 ∫ 𝑓(𝑥)𝑑𝑥と書き直される. これより, 形式的には∞ 𝐵 (ℎ)∞ℎ𝐷𝑘 𝐷𝑘 𝑓(𝑥)𝑑𝑥 = 𝑓(𝑥)𝑑𝑥 + 𝐵𝑘(ℎ) 𝑓 (𝑘−1) (𝑥).
𝐹(𝑥) = 𝑒𝐷𝑒𝑓(𝑥)𝑑𝑥=𝐷 −1 ∫∑ 𝑘! ∫∑ 𝑘!
∫𝑘=0

これより, 整数𝑎 < 𝑏 について, 形式的には
𝑏−1∞ 𝐵 (ℎ)
𝑏
𝑘
𝑓(𝑗
+
ℎ)
=
𝐹(𝑏)
−
𝐹(𝑏)
=
𝑓(𝑥)𝑑𝑥+∑∑ 𝑘! (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)).
∫𝑎
𝑗=𝑎

これは

𝑘=1

𝑘=1

ℎ = 0 とおけば形式的にEuler-Maclaurinの和公式
𝑏−1
∞ 𝐵
𝑏
𝑘
𝑓(𝑗)
=
𝐹(𝑏)
−
𝐹(𝑏)
=
𝑓(𝑥)𝑑𝑥+∑∑ 𝑘! (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)).
∫𝑎
𝑗=𝑎

𝑘=1

に一致する.

2.4 べき乗和(再)𝑓 (𝑚+1) (𝑥) = 0 ならばEuler-Maclaurinの和公式における 𝑅𝑚 が消えて, 𝑓 (𝑚) (𝑥) が定数函数になるので,
𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑚+1 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎))𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∑ 𝑘!
∫𝑎
2
𝑗=𝑎

となる.

𝑘=2

𝑚 を正の整数であるとし, 𝑓(𝑥) = 𝑥𝑚 , 𝑎 = 0, 𝑏 = 𝑛 とすると,(𝑥𝑚 )(𝑘−1) = 𝑚(𝑚 − 1)⋯(𝑚 − 𝑘 + 2)𝑥𝑚−𝑘+1 = 𝑚 𝑘!+ 1 (𝑚 𝑘+ 1)𝑥𝑚+1−𝑘なので,𝐵0 (1) = 𝐵0 = 1, 𝐵1 (1) = 𝐵1 + 1 = 12 , 𝐵𝑛 (1) = 𝐵𝑛 + 𝛿𝑛,1 も使うと,𝑛𝑚𝑚 = 𝑛 𝑥𝑚 𝑑𝑥 + 𝑛𝑚 + 𝐵𝑘 𝑘! 𝑚 + 1 𝑥𝑚+1−𝑘
𝑗
∑
∫0
2 ∑
𝑗=1
𝑘=2 𝑘! 𝑚 + 1 ( 𝑘 )𝑛𝑚 + 1 𝑚 𝑚 + 1 𝐵𝑘𝑛𝑚+1−𝑘= 𝑚𝑛𝑚+1++1 2 𝑚+1 ∑
𝑘=2 ( 𝑘 )𝑚
= 𝑚 +1 1 ∑ (𝑚 𝑘+ 1)𝐵𝑘(1)𝑛𝑚+1−𝑘
𝑘=0
(𝑥)
𝐵
= 𝑚+1𝑚 +−1𝐵𝑚+1 .

最後の等号で

𝑚+1
𝐵𝑚+1(𝑥 + ℎ) = ∑ (𝑚 𝑘+ 1)𝐵𝑘(ℎ)𝑥𝑚+1−𝑘を使った. 以上の結果は前節で得た結果と一致する.

□

𝑘=0

2.5 Stirlingの近似公式

log𝑛! = log1 + log2 + ⋯ + log𝑛 の log𝑥 の積分による近似をEuler-Maclaurinの和公式を使って精密化すると, Stirlingの近似公式

log𝑛! = 𝑛log𝑛 − 𝑛 + 12 log𝑛 + log √⎯2𝜋⎯⎯⎯ + 𝑜(1)の精密化が得られることを説明しよう.

𝑓(𝑥) = log𝑥 にStirlingの和公式を適用すると, 𝑘 = 1,2,3,… に対して,
𝑘−1
𝑓 (𝑘) (𝑥) = (−1) 𝑥(𝑘𝑘 − 1)! .
𝐾 は3以上の奇数とする. 3以上の奇数 𝑘 について 𝐵𝑘 = 0 なので, 𝑁 > 𝑛 のとき,
𝑁

log𝑛! = log𝑁! + log𝑛 − ∑ log𝑗
𝑗=𝑛

𝐾−1
𝑁
1 + 𝑅𝐾,𝑁= log𝑁! + log𝑛 − (∫𝑛 log𝑥𝑑𝑥 + log𝑛 +2 log𝑁 + ∑ 𝑘(𝑘𝐵−𝑘 1) ( 𝑁1𝑘−1 − 𝑛𝑘−1
)
)
𝑘=2
𝐾−1
= log𝑁! − (𝑁 log𝑁 − 𝑁 + 12 log𝑁 ) − ∑ 𝑘(𝑘𝐵−𝑘 1) 𝑁1𝑘−1
𝑘=2
𝐾−1
1 + 𝑅𝐾,𝑁 ,+ 𝑛log𝑛 − 𝑛 + 12 log𝑛 + ∑ 𝑘(𝑘𝐵−𝑘 1) 𝑛𝑘−1
𝑘=2
𝑁 𝐵𝐾̃ (𝑥) (−1)𝐾−1
𝑅𝐾,𝑁 = (−1)𝐾−1 ∫𝑛 𝐾 𝑥𝐾 𝑑𝑥
ただし, 𝐵𝑛̃ (𝑥) = 𝐵𝑛 (⌊𝑥⌋) とおいた.
ここでは, 𝑁 → ∞ のとき
log𝑁! − (𝑁 log𝑁 − 𝑁 + 12 log𝑁 ) → √⎯2𝜋⎯⎯⎯となることは既知であるものとする. 例えば, ノート「10 Gauss積分, ガンマ函数, ベータ函数(http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/10%20Gauss%2C%20Gamma%2C%20Beta.ipynb)」「12Fourier解析 (http://nbviewer.jupyter.org/github/genkuroki/Calculus/blob/master/12%20Fourier%20analysis.ipynb)」のStirlingの近似公式の節を参照して欲しい. 以下ではそれらのノートよりも精密な結果を得る.
このとき, 上の結果で𝑁 → ∞ とすると,

𝐾−1
1 + 𝑅𝐾 ,log𝑛! = 𝑛log𝑛 − 𝑛 + 12 log𝑛 + log √⎯2𝜋⎯⎯⎯ + ∑ 𝑘(𝑘𝐵−𝑘 1) 𝑛𝑘−1
𝑘=2
∞ 𝐵̃ (𝑥) (−1)𝐾−1
1 .
𝑅𝐾 = (−1)𝐾−1 ∫𝑛 𝐾𝐾 𝑥𝐾 𝑑𝑥 = 𝑂 ( 𝑛𝐾−1)
𝐾 = 2𝐿 + 1 とおくことによって次が得られる: 正の整数 𝐿 に対して,
𝐿
1 +𝑂 1 .
log𝑛! = 𝑛log𝑛 − 𝑛 + 12 log𝑛 + log √⎯2𝜋⎯⎯⎯ + ∑ (2𝑙)(2𝑙𝐵2𝑙− 1) 𝑛2𝑙−1( 𝑛2𝐿 )
𝑙=1

これが求めていた結果である.
例えば,In [18]:

Out[18]:

𝐿 = 2 のとき, 𝐵2 = 16 , 𝐵4 = − 301 なので,
1 +𝑂 1 .
log𝑛! = 𝑛log𝑛 − 𝑛 + 12 log𝑛 + log √⎯2𝜋⎯⎯⎯ + 12𝑛1 − 360𝑛3
( 𝑛4 )

1 x = symbols("x")
2 series(exp(x/12-x^3/360), x, n=4)𝑥2 − 139𝑥3 + 𝑂 (𝑥4 )1 + 12𝑥 + 288
51840

2.6 Poissonの和公式とEuler-Maclaurinの和公式の関係
Poissonの和公式とは, 急減少函数

𝑓(𝑥) に対して,̂∑ 𝑓(𝑚) = ∑ 𝑓 (𝑛),
𝑚∈ℤ

𝑛∈ℤ

̂ = ∫ 𝑓(𝑥)𝑒2𝜋𝑖𝑝𝑥 𝑑𝑥
𝑓 (𝑝)ℝが成立するという結果であった. これの右辺は以下のように変形できる:∞̂ = ∑ ∫ 𝑓(𝑥)𝑒2𝜋𝑖𝑛𝑥 𝑑𝑥 = ∫ 𝑓(𝑥)𝑑𝑥 + 2 ∑ ∫ 𝑓(𝑥)cos(2𝜋𝑛𝑥)𝑑𝑥
𝑓 (𝑛)
∑
ℝ
𝑛∈ℤ
𝑛∈ℤ ℝ
𝑛=1 ℝ
∞
= ∫ℝ 𝑓(𝑥)𝑑𝑥 − ∑ ∫ℝ 𝑓 ′ (𝑥) sin(2𝜋𝑛𝑥)
𝜋𝑛 𝑑𝑥
𝑛=1
∞
= ∫ℝ 𝑓(𝑥)𝑑𝑥 + ∫ℝ (− ∑ sin(2𝜋𝑛𝑥)𝑓 ′ (𝑥)𝑑𝑥𝜋𝑛)𝑛=1
2つ目の等号では 𝑒2𝜋𝑖𝑛𝑥 + 𝑒−2𝜋𝑖𝑛𝑥 = 2cos(2𝜋𝑛𝑥) を用い, 3つ目の等号では部分積分を実行し, 4つ目の等号では無限和と積分の順序を交換した. それらの操作は 𝑓(𝑥) が急減少函数であれば容易に正当化される.
一方, Euler-Maclaurinの和公式のを使う場合から,𝐵1 (𝑥 − ⌊𝑥⌋) = 𝑥 − ⌊𝑥⌋ − 12
𝑓(𝑚) = ∫ℝ 𝑓(𝑥)𝑑𝑥 + ∫ℝ (𝑥 − ⌊𝑥⌋ − 12 ) 𝑓 ′ (𝑥)𝑑𝑥∑𝑚∈ℤが導かれる. これは部分積分によって得られる次の公式からただちに導かれる易しい公式であることにも注意せよ:

1 𝑓 ′ (𝑥)𝑑𝑥 = 𝑥 − 𝑛 − 1 𝑓(𝑥) 𝑛+1 − 𝑛+1 𝑓(𝑥)𝑑𝑥 − 𝑛
𝑥
−
𝑛
−
∫𝑛 (
[(
2)
2 ) ]𝑛 ∫𝑛
𝑛+1= 𝑓(𝑛 + 1)2 − 𝑓(𝑛) − ∫𝑛 𝑓(𝑥)𝑑𝑥.
以上の2つの結果を比較すると, Poissonの和公式とEuler-Maclaurinの和公式の 𝐵1 (𝑥 − ⌊𝑥⌋) を使った場合は,∞
𝑥 − ⌊𝑥⌋ − 12 = − ∑ sin(2𝜋𝑛𝑥)
𝜋𝑛
𝑛+1

𝑛=1

という公式で結び付いていることがわかる. この公式を認めれば, Euler-Maclaurinの和公式のPoissonの和公式が導かれる.

(∗)

𝐵1 (𝑥 − ⌊𝑥⌋) を使った場合から∗∗

公式( )の左辺はいわゆるのこぎり波であり, 右辺はそのFourier級数である. 公式( )はFourier級数論における非常に有名な公式
であり, 本質的にそれと同じ公式はFourier級数論について書かれた文献には例として必ず載っていると言ってよいくらいである. (Fourier級数論より, 公式( )は が整数でないときには実際に成立していることがわかる.)∗ 𝑥

このように, のこぎり波のFourier級数展開という非常に特殊な公式はPoissonの和公式という一般的な公式を導くだけの力を持っているのである.
まとめ: のこぎり波のFourier級数展開は部分積分を通してPoissonの和公式と本質的に同値である!

□

この節で解説したことは次の文献で指摘されている:
Tim Jameson, An elementary derivation of the Poisson summation formula(http://www.maths.lancs.ac.uk/jameson/poisson.pdf)
In [19]:



1
2
3
4
5
6
7
8

B_1(x) = x - 1/2
b(x) = B_1(x - floor(x))
S(N,x) = -sum(n->sin(2π*n*x)/(π*n), 1:N)
x = -2:0.001:1.999
N = 10plot(size=(400,200), ylim=(-0.6,1.2), legend=:top)
plot!(x, b.(x), label="B_1(x-[x]) = x - [x] -1/2")
plot!(x, S.(N,x), label="partial sum of Fourier series (N=$N)")

Out[19]:𝐵𝑘(𝑥 − ⌊𝑥⌋) のFourier級数展開の節を見ればわかるように,∞ cos(2𝜋𝑛𝑥) ∞ sin(2𝜋𝑛𝑥), ∑ 𝑛𝑘∑𝑛=1 𝑛𝑘𝑛=1の型のFourier級数の収束先は平行移動と定数倍の違いを除いて周期的Bernoulli多項式になる. □補足: このノートの上の方の周期的Bernoulli多項式

2.7 台形公式とPoissonの和公式の関係
簡単のため

𝑓(𝑥) は ℝ 上の急減少函数であるとし, 𝑎,𝑏 ∈ ℤ かつ 𝑎 < 𝑏 であると仮定する. このとき,
sin(2𝜋𝑛𝑥) ′ = 2cos(2𝜋𝑛𝑥)( 𝜋𝑛 )を使った部分積分とノコギリ波のFourier級数展開∞
1 = 𝐵1 (𝑥 − ⌊𝑥⌋)
− ∑ sin(2𝜋𝑛𝑥)
=
𝑥
−
⌊𝑥⌋
−
𝜋𝑛
2
𝑛=1

より,

𝑏

∞

𝑏

𝑏

𝑓(𝑥)𝑒−2𝜋𝑖𝑛𝑥 𝑑𝑥 = ∫𝑎 𝑓(𝑥)𝑑𝑥 + ∑ 2 ∫𝑎 𝑓(𝑥)cos(2𝜋𝑛𝑥)𝑑𝑥
∑
∫
𝑎
𝑛∈ℤ
𝑛=1
∞ 𝑏
𝑏
= ∫𝑎 𝑓(𝑥)𝑑𝑥 − ∑ ∫𝑎 𝑓 ′ (𝑥) sin(2𝜋𝑛𝑥)
𝜋𝑛 𝑑𝑥
𝑛=1
𝑏
𝑏
= ∫𝑎 𝑓(𝑥)𝑑𝑥 + ∫𝑎 𝑓 ′ (𝑥) (𝑥 − ⌊𝑥⌋ − 12 ) 𝑑𝑥
𝑏−1 𝑚+1𝑏= ∫𝑎 𝑓(𝑥)𝑑𝑥 + ∑ ∫𝑚 𝑓 ′ (𝑥) (𝑥 − ⌊𝑥⌋ − 12 ) 𝑑𝑥
𝑚=𝑎
𝑏−1 1
𝑏
= ∫𝑎 𝑓(𝑥)𝑑𝑥 + ∑ ∫0 𝑓 ′ (𝑡 + 𝑚) (𝑡 − 12 ) 𝑑𝑡
𝑚=𝑎
𝑏−1
1
𝑏
1
= ∫𝑎 𝑓(𝑥)𝑑𝑥 + ∑ ([𝑓(𝑡 + 𝑚) (𝑡 − 12 )] − ∫0 𝑓(𝑡 + 𝑚)𝑑𝑡)0
𝑚=𝑎𝑏−1 𝑓(𝑚) + 𝑓(𝑚 + 1)
=∑
2
𝑚=𝑎
𝑏
= ∑ 𝑓(𝑚) − 𝑓(𝑎) +2 𝑓(𝑏) .
𝑚=𝑎
要するに, ノコギリ波 𝐵1 (𝑥 − ⌊𝑥⌋) のFourier級数展開の公式を認めて使えば,
𝑏−1𝑏𝑏−2𝜋𝑖𝑛𝑥 𝑑𝑥 = 𝑓(𝑚) + 𝑓(𝑚 + 1) = 𝑓(𝑚) − 𝑓(𝑎) + 𝑓(𝑏)
𝑓(𝑥)
𝑒
∑ ∫𝑎
∑
∑
2
2
𝑛∈ℤ

𝑚=𝑎

𝑚=𝑎

という公式が得られる この公式の右辺は積分の近似計算で使われる台形公式に一致しており 途中の計算では

𝐵1 (𝑥 − ⌊𝑥⌋) を

3 ゼータ函数への応用

𝑠 > 1 のとき(より一般には Re𝑠 > 1 のとき),∞
𝜁(𝑠) = ∑ 𝑛1𝑠𝑛=1

は絶対収束しているのであった. これにEuler-Maclaurinの和公式

𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑛 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) + 𝑅𝑛 ,𝑓(𝑗)=𝑓(𝑥)𝑑𝑥
+
∑
∑
∫𝑎
2
𝑗=𝑎
𝑘=2 𝑘!
𝑏 𝐵𝑛 (𝑥 − ⌊𝑥⌋)(𝑛)
𝑅𝑛 = (−1)𝑛−1 ∫𝑎𝑛! 𝑓 (𝑥)𝑑𝑥𝑏を適用してみよう.

3.1 解析接続

Re𝑠 > 1 であるとし, 𝑓(𝑥) = 𝑥−𝑠 とおく. このとき,∞∞∞−𝑠 𝑑𝑥 = 𝑥−𝑠+1 = 𝑎−(𝑠−1) ,𝑓(𝑥)𝑑𝑥=𝑓(𝑏) = 𝑏−𝑠 → 0 (𝑏 → ∞).
𝑥∫𝑎∫1[ −𝑠 + 1 ]1 𝑠 − 1
𝐵𝑘 𝑓 (𝑘−1) (𝑥) = 𝐵𝑘 −𝑠 𝑥−𝑠−𝑘+1 , 𝐵𝑛 (𝑥 − ⌊𝑥⌋) 𝑓 (𝑛) (𝑥) = −𝑠 𝐵𝑛 (𝑥 − ⌊𝑥⌋)𝑥−𝑠−𝑛
(𝑛)
𝑘!
𝑛!
𝑘 (𝑘 − 1)なので, 2以上の整数 𝑛 について,𝑛𝜁(𝑠) = 𝑠 −1 1 + 12 − ∑ 𝐵𝑘𝑘 (𝑘−𝑠− 1) + 𝑅𝑛 ,𝑘=2∞−𝑠
𝑅𝑛 = (−1)𝑛−1 ( 𝑛 ) ∫1 𝐵𝑛 (𝑥 − ⌊𝑥⌋)𝑥−𝑠−𝑛 𝑑𝑥.
積分 𝑅𝑛 は Re𝑠 + 𝑛 > 1 ならば絶対収束している. ゆえに, 複素平面全体に 𝜁(𝑠) を自然に拡張する方法(解析接続する方法)が得られた.

𝑎−1 1
∞ 1
∞ 1
𝑛
=
𝑎
=
𝜁(𝑠)
−
そのものではなく,から始まる無限和𝑠𝑠𝑠 にEuler-Maclaurinの和公式を適用すると,
∑
∑
∑
𝑘=1 𝑛
𝑘=𝑎 𝑛
𝑛=1 𝑛
𝑎−1
1 − 𝑛 𝐵𝑘 −𝑠 + 𝑅𝑛,𝑎 ,𝜁(𝑠) = ∑ 𝑛1𝑠 − 1𝑎1−𝑠+− 𝑠 2𝑎𝑠 ∑𝑛=1
𝑘=2 𝑘𝑎𝑠+𝑘−1 (𝑘 − 1)∞
𝑅𝑛,𝑎 = (−1)𝑛−1 (−𝑠𝑛) ∫𝑎 𝐵𝑛 (𝑥 − ⌊𝑥⌋)𝑥−𝑠−𝑛 𝑑𝑥.
In [20]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

# 上の公式における ζ(s) - R_{n,a} の函数化​
# ζ(s) - R_{n,a} = Σ_{m=1}^{a-1} m^{-s} - a^{1-s}/(1-s) + 1/(2a^s)#- Σ_{k=2}^n B_k/(k a^{s+k-1}) binom(-s,k-1)(k is even)#function ApproxZeta(a, n, s)
ss = float(big(s))
z = zero(ss)z += (a ≤ 1 ? zero(ss) : sum(m->m^(-ss), 1:a-1)) # Σ_{m=1}^{a-1} m^{-s}z += -a^(1-ss)/(1-ss)# -a^{1-s}/(1-s)
n == 0 && return zz += 1/(2*a^ss)# 1/(2a^s)n == 1 && return zz -= sum(k -> BB(k)/(k*a^(ss+k-1))*binom(-ss,k-1), 2:2:n)#-Σ_{k=2}^n B_k/(k a^{s+k-1}) binom(-s,k-1)(k is even)end​
A = ApproxZeta(40, 80, big"0.5")
Z = zeta(big"0.5")@show A
@show Z;A = -1.460354508809586812889499152515298012467229331012581490542886087825530529474572
Z = -1.460354508809586812889499152515298012467229331012581490542886087825530529474503

Re𝑠 > 0 のとき,は𝑎 → ∞ で 0 に収束するので,

1 − 𝑛 𝐵𝑘 −𝑠 + 𝑅𝑛,𝑎2𝑎𝑠 ∑𝑘=2 𝑘𝑎𝑠+𝑘−1 (𝑘 − 1)𝑎−1
1−𝑠
𝜁(𝑠) = 𝑎→∞lim (∑ 𝑛1𝑠 − 1𝑎− 𝑠 ) (Re𝑠 > 0)

𝑛=1
𝑎−1 1
が成立することがわかる. これは, Dirichlet級数の部分和
𝑠 から補正項
∑
𝑛=1 𝑛
𝑎1−𝑠

1−𝑠
を引き去ってから, Dirichlet級数の総和を取れば, 0 < Re𝑠 < 1 でも収束して, 𝜁(𝑠) の正確な値が得られることを意味している.

In [21]:



1
2
3
4
5
6
7
8
9
10

# 上の結果のプロット​
ApproxZeta0(a, s) = sum(n->n^(-s), 1:a-1) - a^(1-s)/(1-s)
a = 100
s = 0.05:0.01:0.95@time z = zeta.(s)@time w = ApproxZeta0.(a, s)plot(size=(400, 250), legend=:bottomleft, xlabel="s")plot!(s, z, label="zeta(s)", lw=2)
plot!(s, w, label="Euler-Maclaurin sum for n=0, a=$a", lw=2, ls=:dash)0.191314 seconds (143.69 k allocations: 9.440 MiB, 7.58% gc time, 44.93% compilation tim
e)
0.055798 seconds (117.52 k allocations: 8.121 MiB, 98.91% compilation time)Out[21]:

In [22]:



1
2
3
4
5
6
7
8
9
10
11
12
13

# さらに項の数を1つ増やした場合のプロット​# ζ(s) - R_{1,a} = Σ_{n=1}^{a-1} n^{-s} - a^{1-s}/(1-s) + 1/(2a^s)#
ApproxZeta1(a, s) = sum(n->n^(-s), 1:a-1) - a^(1-s)/(1-s) + 1/(2*a^s)​s = -0.95:0.01:0.5
a = 10^3@time z = zeta.(s)@time w = ApproxZeta1.(a,s)
plot(size=(400, 250), legend=:bottomleft, xlabel="s")plot!(s, z, label="zeta(s)", lw=2)
plot!(s, w, label="Euler-Maclaurin sum for n=1, a=$a", lw=2, ls=:dash)0.000062 seconds (4 allocations: 1.406 KiB)
0.062979 seconds (114.81 k allocations: 7.860 MiB, 94.18% compilation time)Out[22]:

In [23]:



1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

# さらに一般の場合のプロット## Euler-Maclaurinの和公式で ζ(s) の負の s での値をぴったり近似できていることがわかる.
​[(-m, zeta(-m), Float64(ApproxZeta(2, 17, -m))) for m = 0:12] |> display​
n = 10s = -1.5:0.05:0.5
a = 10@time z = zeta.(s)@time w = ApproxZeta.(a, n, s)
P1 = plot(size=(400, 250), legend=:bottomleft, xlabel="s")plot!(s, z, label="zeta(s)", lw=2)
plot!(s, w, label="Euler-Maclaurin sum for a=$a, n=$n", lw=2, ls=:dash)​n = 17s = -16:0.05:-2.0
a = 2@time z = zeta.(s)@time w = ApproxZeta.(a, n, s)
P2 = plot(size=(400, 250), legend=:topright, xlabel="s")plot!(s, z, label="zeta(s)", lw=2)
plot!(s, w, label="Euler-Maclaurin sum for a=$a, n=$n", lw=2, ls=:dash)13-element Vector{Tuple{Int64, Float64, Float64}}:(0, -0.5, -0.5)(-1, -0.08333333333333338, -0.08333333333333333)(-2, -0.0, -1.2954252832641667e-77)(-3, 0.008333333333333344, 0.008333333333333333)(-4, -0.0, -3.454467422037778e-77)(-5, -0.0039682539682539715, -0.003968253968253968)(-6, -0.0, 0.0)(-7, 0.004166666666666669, 0.004166666666666667)(-8, -0.0, 0.0)(-9, -0.00757575757575758, -0.007575757575757576)(-10, -0.0, -4.421718300208356e-75)(-11, 0.0210927960927961, 0.021092796092796094)(-12, -0.0, 0.0)
0.000034 seconds (4 allocations: 592 bytes)
0.095275 seconds (216.85 k allocations: 13.766 MiB, 91.64% compilation time)
0.000105 seconds (4 allocations: 2.562 KiB)
0.053619 seconds (492.29 k allocations: 19.944 MiB, 27.20% gc time)Out[23]:

In [24]:

1 display(P1)上と下のグラフを見ればわかるように, Euler-Maclaurinの和公式によって負の実数でのる. 実はを実部が負の複素数まで拡張してもこの近似はうまく行っている.

𝜁(𝑠)In [25]:𝜁 函数の値を非常によく近似できてい1 display(P2)3.2 ζ(2)の近似計算∞
𝜁(2) = ∑ 𝑛12 を計算せよという問題はBasel問題と呼ばれているらしい. Basel問題はEulerによって1743年ころに解かれたら
𝑛=1

しい. Eulerがどのように考えたかについては次の文献を参照せよ.

杉本敏夫, バーゼル問題とオイラー (http://www.kurims.kyoto-u.ac.jp/~kyodo/kokyuroku/contents/pdf/1583-12.pdf), 2007年8月23日, 数理解析研究所講究録, 第1583巻, 2008年, pp.159-167
Eulerは

𝜁(2) の近似値を自ら開発したEuler-Maclaurinの和公式を使って精密に計算したらしい.

近似式

𝑎−1
𝑛
1−𝑠
−𝑠
𝑘
𝜁(𝑠) ≈ ∑ 𝑛1𝑠 − 1𝑎− 𝑠 + 2𝑎1𝑠 − ∑ 𝑘𝑎𝐵𝑠+𝑘−1(𝑘 − 1)
𝑛=1

を用いて,
る.

𝑘=2

𝜁(2) を計算してみよう. 3以上の奇数 𝑛 について 𝐵𝑛 = 0 となるので, 𝑛 = 2𝑚 のとき, 右辺の項数は 𝑎 + 𝑚 + 1 にな

𝑎 = 10, 𝑚 = 9 とし, 20項の和を取ると,
𝜁(2) ≈ 1.64493 40668 4749⋯
𝜋2 = 1.64493 40668 4822⋯ と小数点以下第11桁まで一致している.
となり, 正確な値
6𝜋2 を得る. Eulerは競争相手に議論に厳密性に欠けるとして様々な批判を受けたのだが, 以上のような数値
Eulerは後に 𝜁(2) =6
計算の結果を知っていたので, 正解を得たという確信は微塵も揺らがなかっただろうと思われる.
例えば,

注意: 論理的に厳密な証明の方法が発達した現代においても, 人間は常に証明を間違う可能性がある. 人間が行った証明は絶対的には信用できない. だから, たとえ証明が完成したと思っていたとしても, 可能ならば数値計算によって論理的に厳密な証明以外の証拠を作っていた方が安全だと思われる.

□

注意: 数学のノートを作りながら, 気軽に数値的証拠も同時に得るための道具として, 筆者がこのノート作成のために用いている
Julia言語 (https://julialang.org/)とJupyter (http://jupyter.org/)とNbextensions (https://github.com/ipythoncontrib/jupyter_contrib_nbextensions)のLive Markdown Previewはこれを書いている時点で相当に優秀な道具であるように思われる.

□

In [26]:



1
2
3
4
5
6
7
8
9
10
11
12

# 20項の和​N = 20[(m, N-m-1, 2m, ApproxZeta(N-m-1, 2m, 2) - big(π)^2/6) for m in 2:N÷2-1] |> display​
m = 9
a = N-m-1
Z = big(π)^2/6
A = ApproxZeta(a, m, 2)@show a,m@show Z
@show A;

8-element Vector{Tuple{Int64, Int64, Int64, BigFloat}}:(2, 17, 4, -5.7745179386347483379778894047835850369958540757839935768100161932366414526175
8e-11)(3, 16, 6, 4.80812735239562509501346011215032587838986605495815313740843048777477650932482
9e-13)(4, 15, 8, -8.6308875139430442246152364654652069706509110461365277080262921867238658169664
92e-15)(5, 14, 10, 3.1162179785273853280542355730238711734665867971863964368976627204145528522974
51e-16)(6, 13, 12, -2.200847274100542514575619216657515798396053843860275532661691594239630209744
406e-17)(7, 12, 14, 3.0352489438578151477776773837113166940196569351033194321813552488718204060625
84e-18)(8, 11, 16, -8.335321043122531064769674746337938450627967961329547742403411230422546148897
753e-19)(9, 10, 18, 4.7466018143920053127140275780279703065395409350513421647372241615147960630210
67e-19)(a, m) = (10, 9)
Z = 1.644934066848226436472415166646025189218949901206798437735558229370007470403185
A = 1.644934066847493071302595112118921642731166540690350214159737969261778785588307

3.3 s = 1でのζ(s)の定数項がEuler定数になること∞
𝜁(𝑠) = ∑ 𝑛1𝑠 にEuler-Maclaurinの和公式を使って, 2以上の 𝑛 について次の公式が得られるのであった:
𝑛=1

𝑛 = 1 の場合には𝑛𝜁(𝑠) = 𝑠 −1 1 + 12 − ∑ 𝐵𝑘𝑘 (𝑘−𝑠− 1) + 𝑅𝑛 ,𝑘=2∞−𝑠
𝑅𝑛 = (−1)𝑛−1 ( 𝑛 ) ∫1 𝐵𝑛 (𝑥 − ⌊𝑥⌋)𝑥−𝑠−𝑛 𝑑𝑥.

𝑏

𝑏

𝑏

𝑏

𝑏−1

∑ 𝑓(𝑗) = ∫𝑎 𝑓(𝑥)𝑑𝑥 + 𝑓(𝑎) + ∫𝑎 (𝑥 − ⌊𝑥⌋)𝑓 ′ (𝑥)𝑑𝑥
𝑗=𝑎1= ∫𝑎 𝑓(𝑥)𝑑𝑥 + 𝑓(𝑎) + ∑ ∫0 𝑥𝑓 ′ (𝑥 + 𝑗)𝑑𝑥
𝑗=𝑎を

𝑓(𝑥) = 𝑥−𝑠 , 𝑓 ′ (𝑥) = −𝑠𝑥−𝑠−1 , 𝑎 = 1, 𝑏 = ∞ の場合に適用して,∞ 1
𝜁(𝑠) = 𝑠 −1 1 + 1 − 𝑠 ∑ ∫0 (𝑥 +𝑥𝑗)𝑠+1 𝑑𝑥
𝑗=1を得る. したがって,

1 = 1 − ∞ 1 𝑥 𝑑𝑥.
𝜁(𝑠)
−
lim
𝑠→1 (
𝑠 − 1 ) ∑ ∫0 (𝑥 + 𝑗)2
𝑗=1

𝑥 = 𝑡 − 𝑗 と置換すると,
𝑗+1
𝑗+1
1
− ∫0 (𝑥 +𝑥 𝑗)2 𝑑𝑥 = − ∫𝑗 −(𝑡𝑡2− 𝑗) 𝑑𝑡 = −[log𝑡 + 𝑗𝑡 ]𝑗𝑗= −log(𝑗 + 1) + log𝑗 − 𝑗 + 1 + 1 = 𝑗 +1 1 + log𝑗 − log(𝑗 + 1)なので, これを 𝑗 = 1 から 𝑗 = 𝑁 − 1 まで足し上げることによって,
そして,

𝑛−1 1𝑁
1 − ∑ ∫0 (𝑥 +𝑥 𝑗)2 𝑑𝑥 = ∑ 1𝑗 − log𝑁.
𝑗=1

これの

𝑗=1

𝑁 → ∞ での極限はEuler定数 𝛾 = 0.5772⋯ の定義であった. 以上によって次が示された:
1lim
𝑠→1 (𝜁(𝑠) − 𝑠 − 1 ) = 𝛾 = 0.5772⋯.

3.4 負の整数におけるゼータ函数の特殊値の計算

3 以上の整数 𝑘 について 𝐵𝑘 = 0 なので, 以下の公式で 𝑘 は偶数のみを動くとしてよい:𝑏𝑏𝑓(𝑎) + 𝑓(𝑏) + 𝑚 𝐵𝑘 (𝑓 (𝑘−1) (𝑏) − 𝑓 (𝑘−1) (𝑎)) + 𝑅𝑚 ,𝑓(𝑛)=𝑓(𝑥)𝑑𝑥
+
∑
∑
∫𝑎
2
𝑛=𝑎
𝑘=2 𝑘!
𝑏 𝐵̃ (𝑥)𝑚 𝑓 (𝑚) (𝑥)𝑑𝑥.
𝑅𝑛 = (−1)𝑚−1 ∫𝑎 𝑚!
̃ (𝑥) = 𝐵𝑚(𝑥 − ⌊𝑥⌋) とおいた.
ここで 𝐵𝑚Euler-Maclaurinの和公式を 𝑓(𝑥) = 𝑛−𝑠 , 𝑎 = 1, 𝑏 = ∞ の場合に適用することによって 𝜁(𝑠) は次の形で ℜ𝑠 > 1 − 𝑚 まで自Euler-Maclaurinの和公式:

然に延長(解析接続)されるのであった:𝑚𝑏
𝜁(𝑠) = 𝑠 −1 1 + 12 − 1 −1 𝑠 ∑ (1 −𝑘 𝑠)𝐵𝑘 + (−1)𝑚−1 ∫𝑎 (−𝑠𝑚)𝐵𝑚̃ (𝑥)𝑥−𝑠−𝑚 𝑑𝑥.
𝑘=2

この公式と

𝑘 ≧ 2 のとき (𝑘1) = 0 となることより,

𝜁(0) = 0 −1 1 + 12 = − 12 .
𝑟 は正の整数であるとする. このとき, 𝑚 > 𝑟 とすると (𝑚𝑟 ) = 0 となるので, 𝐵0 = 1, 𝐵1 = −1/2 なので,
𝑟+1
𝜁(−𝑟) = − 𝑟 +1 1 + 12 − 𝑟 +1 1 ∑ (𝑚 𝑘+ 1)𝐵𝑘
𝑘=2
𝑟+1
= − 𝑟 +1 1 ∑ (𝑚 𝑘+ 1)𝐵𝑘 = − 𝑟𝐵+𝑟+11 .
𝑘=0
𝑟 𝑟+1
最後の等号で, Bernoulli数を帰納的に計算するために使える公式
𝐵𝑘 = 0 を用いた. 例えば, 𝑟 = 1 のとき∑
𝑘=0 ( 𝑘 )𝐵0 + 2𝐵1 = 1 + 2(−1/2) = 0 となり, 𝑟 = 2 のとき, 𝐵0 + 3𝐵1 + 3𝐵2 = 1 + 3(−1/2) + 3(1/6) = 0 となる.
以上によって次が証明された:

𝜁(0) = − 12 , 𝜁(−𝑟) = − 𝑟𝐵+𝑟+11 (𝑟 = 1,2,3,…).
これらの公式は 𝐵𝑛 (1) = 𝐵𝑛 + 𝛿𝑛,1 , 𝐵1 = −1/2 を使うと,
𝜁(−𝑟) = − 𝐵𝑟𝑟+1+(1)1 (𝑟 = 0,1,2,…)の形にまとめられる.

3.5 発散級数の有限部分と ζ(-r) の関係
前節の結果

𝜁(−𝑟) = − 𝐵𝑟𝑟+1+(1)1 (𝑟 = 0,1,2,…) は
1 + 1 + 1 + 1 + ⋯ = − 12 ,1 + 2 + 3 + 4 + ⋯ = − 121のような印象的な形式で書かれることもある. ただし, その場合には左辺が通常の無限和ではなく, ゼータ函数の意味であることを了解しておかなければいけない.

𝜁(𝑠) の解析接続

実はさらに解析接続として理解するだけではなく, 「左辺の発散する無限和から適切に無限大を引き去れば右辺に等しくなる」というようなタイプの命題をうまく作ることもできる. 以下ではそのことを解説しよう.

𝜂ℝℝ 上の急減少函数とは ℝ 上の 𝐶 ∞ 函数でそれ自身およ|𝑥| → ∞ で 0 に収束するもののことである.) さらに,
𝜂(0) = 1, 𝜂′(0) = 0

以下, は非負の実数に値を持つ上の急減少函数であると仮定する. (びそのすべての階数の導函数に任意の多項式函数をかけたものが

𝜂(𝑥) = 𝑒−𝑥2 はそのような函数の例になっている.
このとき, 𝜂(𝑥) が急減少函数であることより, 𝑁 > 0 のとき, 級数
∞
𝑟
𝑟
𝑟
𝑟
∑ 𝑛 𝜂(𝑛/𝑁) = 1 𝜂(1/𝑁) + 2 𝜂(2/𝑁) + 3 𝜂(3/𝑁) + ⋯と仮定する. 例えば

𝑟

𝑛=1

𝑁→∞
𝑁→∞

1𝑟 + 2𝑟 + 3𝑟 + ⋯𝑟+1
𝐶𝑁 𝐶 𝜂 𝑟

は常に絶対収束する. が非負の整数のとき,とすると, この級数は発散級数になってしまう. 以下の目標は, Euler-Maclaurinの和公式を使うと, そのでの発散部分が( は と で具体的に決まる定数) の形にまとまることを示すことである. そして, 残った有限部分は常にに収束することも示される.

𝜁(−𝑟)𝐵𝑛̃ (𝑥) = 𝐵𝑛 (𝑥 − ⌊𝑥⌋) と書くことにする.
このとき, 𝑓(𝑥) = 𝜂(𝑥/𝑁) にEuler-Maclaurinの和公式を適用すると, 𝑓(0) = 1, 𝑓 ′ (0) = 𝑓(∞) = 𝑓 ′ (∞) = 0 より,∞∞
1 + ∑ 𝜂(𝑥/𝑁) = ∑ 𝜂(𝑥/𝑁)
𝑛=1
𝑛=0
∞
∞
= ∫0 𝜂(𝑥/𝑁)𝑑𝑥 + 12 + 𝐵2 (𝑓 ′ (∞) − 𝑓 ′ (0)) − ∫0 𝐵2̃2!(𝑥) 𝑁12 𝜂″(𝑥/𝑁)𝑑𝑥∞∞″= 𝑁 ∫0 𝜂(𝑦)𝑑𝑦 + 12 − 𝑁1 ∫0 𝐵2̃ (𝑁𝑦)
2! 𝜂 (𝑦)𝑑𝑦.
ゆえに, 𝜁(0) = −1/2 を使うと,∞∞
𝜂(𝑥/𝑁)−𝑁∑∫0 𝜂(𝑦)𝑑𝑦 = 𝜁(0) + 𝑂(1/𝑁).
𝑛=1

∞

∞

𝑁 → ∞ で発散級数 1 + 1 + 1 + 1 + ⋯ になる無限和 ∑ 𝜂(𝑥/𝑁) から, その発散部分 𝑁 ∫0 𝜂(𝑦)𝑑𝑦 を引き去って,
𝑛=1
𝑁 → ∞ の極限を取ると, 𝜁(0) に収束することを意味している. これが欲しい結果の1つ目である.
𝑟 は正の整数であるとし, 𝑓(𝑥) = 𝑥𝑟 𝜂(𝑥/𝑁) とおく. そのとき, Leibnitz則
𝑟(𝜑(𝑥)𝜓(𝑥))(𝑚) = ∑ (𝑚𝑖)𝜑(𝑖) (𝑥)𝜓 (𝑚−𝑖)(𝑥)
これは

𝑖=0

を使うと,

𝑓 (𝑟+2) (𝑥) = 𝑁12 𝐹(𝑥/𝑁), 𝐹(𝑦) = (𝑟 +0 2)𝑦𝑟 𝜂(𝑟+2) (𝑦) + ⋯ + (𝑟 +𝑟 2)𝑟!𝜂(𝑦)その 𝑓(𝑥) にEuler-Maclaurinの和公式を適用すると, 𝑓 (𝑘) (∞) = 𝑓 (𝑘) (∞) = 0 および,
𝑓(0) = 𝑓 ′ (0) = ⋯ = 𝑓 (𝑟−1) (0) = 𝑓 (𝑟+1) (0) = 0, 𝑓 (𝑟) (0) = 𝑟!
より,∞∞
̃ (𝑥) 𝑓 (𝑟+2) (𝑥)𝑑𝑥
𝐵𝑟+1 𝑟! − 𝐵𝑟+2 0 + (−1)𝑟+1 ∞ 𝐵𝑟+2𝑟 𝜂(𝑛/𝑁) = 𝑓(𝑛) =𝑓(𝑥)𝑑𝑥−𝑛∑∑∫0∫0 (𝑟 + 2)!
(𝑟 + 1)! (𝑟 + 2)!
𝑛=1
𝑛=0
∞
∞ ̃ (𝑁𝑦)= 𝑁 𝑟+1 ∫0 𝑦𝑟 𝜂(𝑦)𝑑𝑦 − 𝑟𝐵+𝑟+11 + (−1)𝑟+1 𝑁1 ∫0 𝐵(𝑟𝑟+2
𝐹(𝑦)𝑑𝑦+2)!
∞= 𝑁 𝑟+1 ∫0 𝑦𝑟 𝜂(𝑦)𝑑𝑦 − 𝑟𝐵+𝑟+11 + 𝑂(1/𝑁).
𝐵𝑟+1
ゆえに, 𝜁(−𝑟) = −𝑟 + 1 を使うと,∞∞∞∑ 𝑛𝑟 𝜂(𝑛/𝑁) − 𝑁 𝑟+1 ∫0 𝑦𝑟 𝜂(𝑦)𝑑𝑦 = 𝜁(−𝑟) + 𝑂(1/𝑁).
𝑛=1

∞

∞

𝑁 → ∞ で発散級数 1𝑟 + 2𝑟 + 3𝑟 + 4𝑟 + ⋯ になる無限和 ∑ 𝑛𝑟 𝜂(𝑥/𝑁) から, その発散部分 𝑁 𝑟+1 ∫0 𝑦𝑟 𝜂(𝑦)𝑑𝑦1′
注意: 以上の計算のポイントは, 非負の急減少函数 𝜂(𝑥) で 𝜂(0) = 1, 𝜂 (0) = 0 を満たすもので発散級数を正則化して得られる級数の場合には, Euler-Maclaurinの和公式の「途中の項」がほとんど消えてしまうことである. 𝐶𝑁 𝑟+1 型の発散項と定数項と𝑂(1/𝑁) の部分の3つの項しか生き残らない. □これは

注意: 以上の結果に関するより進んだ解説については次のリンク先を参照せよ:
Terence Tao, The Euler-Maclaurin formula, Bernoulli numbers, the zeta function, and real-variable analytic continuation(https://terrytao.wordpress.com/2010/04/10/the-euler-maclaurin-formula-bernoulli-numbers-the-zeta-function-and-realvariable-analytic-continuation/), Blog: What's new, 10 April, 2010.
このブログ記事はかなり読み易い.

□

問題: 以上の結果を数値計算でも確認してみよ.

□𝜂(𝑥) = 𝑒−𝑥2 の場合を試してみよ. そのとき,∞∞𝑟 𝜂(𝑦)𝑑𝑦 =𝑟 𝑒−𝑦2 𝑑𝑦 = 1 Γ 𝑟 + 1
𝑦
𝑦
∫0
∫0
2 ( 2 )となっている. □ヒント:

解答例: 次のリンク先のノートを見よ.
黒木玄, ζ(s) の Re s ＜ 1 での様子 (http://nbviewer.jupyter.org/gist/genkuroki/8b13fc9c05bfd669c345b940066b896b)
In [27]:

Out[27]:

1 y = symbols("y", real=true)
2 r = symbols("r", positive=true)
3 integrate(y^r*exp(-y^2), (y, 0, oo))Γ ( 2𝑟 + 12 )2□