学びのTips数学編20 ver.2 2020年5月
円周率、楕円積分、算術幾何平均
円周率 π とは、円周の長さが直径の何倍かを表す数です。π は代数方程式の解にならない超越数 で、とても不思議な数です。π の近似値を求めることは、紀元前 2000 年頃の古代バビロニアに始 まり、これまで様々な方法が提案されてきました。一方、楕円積分は楕円の周の長さを求める際に 現れる積分で、自然科学や工学の様々な分野と関係があります。また、算術幾何平均は相加相乗平 均とも呼ばれ、高校からなじみのあるものです。今回の Tips では、円周率、楕円積分、算術幾何 平均の間に成り立つ美しい公式を紹介します。そして、その公式を用いた、現在知られている円周 率の最速の計算法である算術幾何平均を用いた方法について紹介します。1
楕円の周の長さと完全楕円積分
0 < p < q に対して、 x2 p2 + y2 q2 = 1 で定まる楕円の弧の長さを求めます。第 1 象限にある弧の長さを ℓ とすると、楕円の弧の長さは 4ℓ となるので、ℓ の値を求めればよいです。第 1 象限にある弧は、パラメータ θ を用いて、次のよう に表示できます。 x = p cos θ, y = q sin θ, 0≤ θ ≤ π 2 このとき、弧長を求める公式を用いると、 ℓ = Z π/2 0 s dx dθ 2 + dy dθ 2 dθ = q Z π/2 0 p 1− k2sin2θ dθ ここで、k = s 1− p q 2 としました。0 < p < q より、0 < k < 1 となります。 ここで、0 < k < 1 を満たす実数 k に対して、 E(k) = Z π/2 0 p 1− k2sin2θ dθ とします。E(k) を第 2 種完全楕円積分といいます。また、0 < k < 1 を満たす実数 k に対して、 K(k) = Z π/2 0 1 p 1− k2sin2θ dθ を第 1 種完全楕円積分といいます。定理 1 次が成り立ちます。 dE(k) dk = E(k)− K(k) k , dK(k) dk = E(k) (1− k2)k − K(k) k 証明 dE(k) k = d dk Z π/2 0 p 1− k2sin2θ dθ = Z π/2 0 ∂ ∂k p 1− k2sin2θ dθ = Z π/2 0 −k sin2θ p 1− k2sin2θdθ = 1 k Z π/2 0 1− k2sin2θ p 1− k2sin2θdθ− 1 k Z π/2 0 dθ p 1− k2sin2θ = E(k)− K(k) k dK(k) k = d dk Z π/2 0 dθ p 1− k2sin2θ = Z π/2 0 ∂ ∂k 1 p 1− k2sin2θ ! dθ = Z π/2 0 k sin2θ dθ (1− k2sin2θ)3/2 = Z π/2 0 k 1− k2 d dθ − cos θ p 1− k2sin2θ ! sin θ dθ = k 1− k2 " −pcos θ sin θ 1− k2sin2θ #π/2 0 + k 1− k2 Z π/2 0 cos2θ p 1− k2sin2θdθ = k 1− k2 Z π/2 0 1 k2 1− k2sin2θ− (1 − k2) p 1− k2sin2θ dθ = 1 k(1− k2) Z π/2 0 p 1− k2sin2θ dθ− (1 − k2) Z π/2 0 1 p 1− k2sin2θdθ ! = 1 k(1− k2)(E(k)− (1 − k 2)K(k)) □ a > b > 0 のとき、T (a, b) を T (a, b) = Z π/2 0 dθ √ a2cos2θ + b2sin2θ と定義します。 補題 1 次が成り立ちます。 T (a, b) = T a + b 2 , √ ab 証明 T (a, b) = Z π/2 0 1 √ a2+ b2tan2θ dθ cos θ u = b tan θ と置換します。このとき、 du dθ = b 1 cos2θ = b(1 + tan 2 θ) = b + u 2 b より、 dθ = b b2+ u2du
また、 cos2θ = 1 1 + tan2θ = b2 b2+ u2 より cos θ = √ b b2+ u2 よって、 T (a, b) = Z ∞ 0 du p (a2+ u2)(b2+ u2) 従って、 T a + b 2 , √ ab = 1 2 Z ∞ −∞ du rn a+b 2 2 + u2o(ab + u2) u = 12(v− abv) と置換します。v が 0 から∞ まで動くとき、u は −∞ から ∞ まで単調に増加しま す。また、 du dv = ab + v2 2v2 より、 du = ab + v 2 2v2 dv また、 a + b 2 2 + u2 = (a 2+ v2)(b2+ v2) 4v2 , ab + u 2 = (ab + v2)2 4v2 よって、 T a + b 2 , √ ab = Z ∞ 0 dv p (a2+ v2)(b2+ v2) よって、補題 1 の式が成り立ちます。 □ 補題 2 a > b > 0 のとき、次が成り立ちます。 T (a, b) = 1 aK s 1− b a 2 問 1 補題 2 が成り立つことを確認してみよう。 補題 3 0 < k < 1 のとき、次が成り立ちます。 K(k) = 2 1 +√1− k2K 1−√1− k2 1 +√1− k2
証明 補題 2 において、a = 1, b =√1− k2と置くと、 T (1,√1− k2) = K(k) となります。補題 1 より、 T (1,√1− k2) = T 1 +√1− k2 2 , (1− k 2 )1/4 補題 2 より、 T 1 +√1− k2 2 , (1− k 2)1/4 = 2 1 +√1− k2K 1−√1− k2 1 +√1− k2 よって、補題 3 の式が成り立ちます。 □ 補題 4 0 < k < 1 のとき、次が成り立ちます。 E(k) = 1 +√1− k2E 1−√1− k2 1 +√1− k2 −√1− k2K(k) 証明 k′ =√1− k2, ℓ = 1− k′ 1 + k′ とおきます。k ′, ℓ を k について微分すると、 dk′ dk =− k k′ dℓ dk = dℓ dk′ dk′ dk = 2k2 (1 + k′)2k′k = 2ℓ kk′ 補題 3 より、 (1 + k′)K(k) = 2K(ℓ) 両辺を k で微分すると、 (1 + k′) d dkK(k) + dk′ dkK(k) = 2 d dℓK(ℓ) dℓ dk (1) 定理 1 より、(1) の左辺は、 1 + k′ kk′2 (E(k)− k ′2K(k))− k k′K(k) (2) 定理 1 より、(1) の右辺は、 4 E(ℓ) (1− ℓ2)kk′ − 4 K(ℓ) kk′ ここで、 1− ℓ2 = 1− (1− k ′)2 (1 + k′)2 = 4k′ (1 + k′)2 よって、 4 E(ℓ) (1− ℓ2)kk′ = (1 + k′)2 k′2k E(ℓ)
補題 3 より、 4 kk′K(ℓ) = 2 1 + k′ kk′ K(k) よって、(1) の右辺は、 (1 + k′)2 kk′2 E(ℓ)− 2(1 + k′) kk′ K(k) (3) (2) と (3) より、E(k) について整理すると、補題 4 の式を得ます。 □ a > b > 0 のとき、 E(a, b) = Z π/2 0 p a2cos2θ + b2sin2θ dθ とおきます。 補題 5 a > b > 0 のとき、次が成り立ちます。 E(a, b) = aE s 1− b a 2 問 2 補題 5 が成り立つことを確認してみよう。 補題 6 a > b > 0 のとき、次が成り立ちます。 2E a + b 2 , √ ab
− E(a, b) = abT (a, b)
証明 具体的に計算すると、 v u u t1− 2√ab a + b !2 = a− b a + b となることが分かります。よって、補題 5 より、 E a + b 2 , √ ab = a + b 2 E a− b a + b (4) となります。k = s 1− b a 2 とおきます。このとき、√1− k2 = b/a となるので、 1−√1− k2 1 +√1− k2 = a− b a + b が成り立ちます。よって、 a + b 2 E a− b a + b = a + b 2 E 1−√1− k2 1 +√1− k2 (5)
(4) と (5) より、 E a + b 2 , √ ab = a + b 2 E 1−√1− k2 1 +√1− k2 一方、補題 5 より、E(a, b) = aE(k) となります。よって、 2E a + b 2 , √ ab − E(a, b) = (a + b)E 1−√1− k2 1 +√1− k2 − aE(k) (6) 補題 4 より、 E 1−√1− k2 1 +√1− k2 = E(k) + √ 1− k2K(k) 1 +√1− k2 (7) √ 1− k2 = b/a であるから、 E(k) +√1− k2K(k) 1 +√1− k2 = aE(k) + bK(k) a + b (8) (6) と (7) と (8) より、 2E a + b 2 , √ ab
− E(a, b) = bK(k) = abT (a, b)
最後の等式は、補題 2 を用いました。 □
2
算術幾何平均
0≤ b ≤ a を満たす実数 a, b に対して、数列 {an}, {bn} を次で定義します。 a0 = a, b0 = b, an+1= an+ bn 2 , bn+1 = p anbn, (n = 0, 1, 2,· · · ) (9) 補題 7 任意の n≥ 0 に対して、bn ≤ anが成り立ちます。 証明 数学的帰納法で示します。n = 0 のときは、定義より b0 ≤ a0 となります。n = k のとき bk≤ akが成り立つとします。 ak+1− bk+1 = ak+ bk 2 − p akbk = (√ak− √ bk)2 2 ≥ 0 となるので、n = k + 1 のときも主張が成り立ちます。よって、任意の n≥ 0 に対して、bn≤ anが 成り立ちます。 □ 補題 8 任意の n≥ 0 に対して、an+1 ≤ an, bn≤ bn+1が成り立ちます。 証明 補題 7 より、 an− an+1 = an− an+ bn 2 = an− bn 2 ≥ 0 bn+1− bn= p anbn− bn= p bn( √ an− p bn)≥ 0となり、題意は成り立ちます。 □ 補題 7, 8 より、 b = b0 ≤ b1 ≤ · · · ≤ bn≤ an ≤ · · · ≤ a1 ≤ a0 = a となります。有界な単調数列は収束するので、数列{an}, {bn} は n → ∞ のときに収束します。an の極限値を α、bnの極限値を β とします。β ≤ α となります。cn= p a2 n− b2nとします。 補題 9 任意の n≥ 1 に対して、 0≤ cn ≤ c0 2n 証明 まず、n≥ 1 において、 cn = p a2 n− b2n = s an−1+ bn−1 2 2 − an−1bn−1 = r (an−1− bn−1)2 4 = an−1− bn−1 2 が成り立ちます。上の式で n = 1 とすると、 c1 = a0− b0 2 = a− b 2 ≤ √ a2− b2 2 = c0 2 また、{bn} は単調増加な数列であることを用いると、n ≥ 2 において次の不等式が成り立ちます。 cn= an−1− bn−1 2 = 1 2 an−2+ bn−2 2 − bn−1 ≤ 1 2 an−2+ bn−2 2 − bn−2 = 1 2 an−2− bn−2 2 = cn−1 2 よって、 0≤ cn ≤ cn−1 2 ≤ cn−2 22 ≤ · · · ≤ c0 2n □ 補題 9 より、n → ∞ のとき、cn → 0 となるので、β = α であることが分かります。α は a, b か
ら定まる実数であり、これを AGM(a, b) と書きます。an, bnが AGM(a, b) に収束する速度は非常に
高速であることが知られています。
問 3 AGM(a, 0) = 0, AGM(a, a) = a であることを確かめてみよう。 命題 1 a > 0 のとき、AGM(a, b) = aAGM(1, b/a) が成り立ちます。 証明 a1 = a + b 2 , b1 = √ ab, an+1 = an+ bn 2 , bn+1 = p anbn, n = 1, 2, . . . とします。また、 c1 = 1 + b/a 2 , d1 = r 1· b a, cn+1 = cn+ dn 2 , dn+1= p cndn, n = 1, 2, . . . とします。このとき、任意の自然数 n に対して、an = acn, bn = adnとなることを数学的帰納法で 示します。n = 1 のときは明らかに成り立ちます。n = k のとき成り立つとします。 ak+1 = ak+ bk 2 = ack+ adk 2 = a ck+ dk 2 = ack+1
bk+1 = p akbk = p ackadk = a p ckdk = adk+1 よって、n = k +1 のときも成り立ちます。以上より、任意の自然数 n に対して、an= acn, bn = adn となることが分かりました。よって、 AGM(a, b) = lim
n→∞an= a limn→∞cn= aAGM(1, b/a)
□
3
円周率、楕円積分、算術幾何平均の関係
ガウスは次の定理を示しました。 定理 2 0 < k < 1 のとき、次が成り立ちます。 K(k) = π 2 · 1 AGM(1,√1− k2) 円周率、楕円積分、算術幾何平均が結びついた大変きれいな定理です。この定理を示します。 証明 a = 1, b = √1− k2としたとき、(9) のように a n, bnを定義します。α = AGM(1, √ 1− k2) とします。n→ ∞ のとき、an, bn → α となります。補題 1 と補題 2 より、 K(k) = T (1,√1− k2) = T (a 1, b1) = T (a2, b2) = · · · = T (α, α) = 1 α Z π/2 0 dθ = π 2 · 1 AGM(1,√1− k2) □4
円周率の超高速計算法
算術幾何平均を用いることで、円周率を非常に高速に計算する公式を導くことができます。こ の公式を示すために次の命題を用います。 命題 2 (ルジャンドルの関係式) 任意の 0 < k < 1 に対して、次が成り立ちます。 E(k) K(k)+ E(√1− k2) K(√1− k2) − 1 = π 2 · 1 K(k)K(√1− k2) 証明 k′ =√1− k2とします。 (E(k)− K(k))K(k′) + E(k′)K(k) = π 2 (10) を示します。 d dkk ′ =−k k′が成り立ちます。また、√1− k′2 = k が成り立ちます。定理 1 より、 d dk{(E(k) − K(k))K(k ′) + E(k′)K(k)} = d dk(E(k)− K(k)) K(k′) + (E(k)− K(k)) d dkK(k ′) + d dkE(k ′)K(k) + E(k′) d dkK(k) = E(k)− K(k) k − E(k)− k′2K(k) kk′2 K(k′) + (E(k)− K(k))E(k ′)− k2K(k′) k′k2 −k k′ +E(k ′)− K(k′) k′ −k k′ K(k) + E(k′)E(k)− k ′2K(k) kk′2 = 0 よって、(E(k)− K(k))K(k′) + E(k′)K(k) は k によらず一定です。 K(0) = E(0) = π 2, E(1) = 1 が成り立ちます。0 < k < 1 のとき、E(k)≤ K(k) が成り立ちます。よって、0 < k < 1 のとき、 0 ≤ (K(k) − E(k))K(k′) = Z π/2 0 k2sin2θ p 1− k2sin2θ dθ ! Z π/2 0 dθ √ cos2θ + k2sin2θ ! ≤ k Z π/2 0 1 p 1− k2sin2θ dθ ! Z π/2 0 kdθ p k2(cos2θ + sin2θ) ! ≤ k · K(k) · π 2 k → 0 とすると、k · K(k) ·π2 → 0 となります。よって、 lim k→0(K(k)− E(k))K(k ′) = 0 となります。よって、 lim k→0{(K(k) − E(k))K(k ′) + E(k′)K(k)} = E(1)K(0) = π 2 よって、(10) が成り立ちます。(10) の両辺を K(k)K(k′) で割ると求める式を得ます。 □ 命題 3 0 < k < 1 とします。a = 1, b =√1− k2としたとき、(9) のように a n, bnを定義します。 cn= p a2 n− b2nとします。このとき、次が成り立ちます。 E(k) K(k) = 1− ∞ X n=0 2n−1c2n 証明 補題 1 と補題 6 より、 2n+1(E(an+1, bn+1)− a2n+1T (an+1, bn+1))− 2n(E(an, bn)− a2nT (an, bn)) = 2n(2E(an+1, bn+1)− E(an, bn)) + (2na2n− 2n+1a2n+1)T (an, bn) = 2n−1(2anbn− 4a2n+1+ 2a 2 n)T (an, bn) = 2n−1(a2n− b 2 n)T (an, bn) = 2n−1c2nK(k) 最後の等式は、任意の自然数 n に対して、 T (an, bn) = T (an−1, bn−1) =· · · = T (a, b) = K(k) (11)
が成り立つことを用いました。T (a, b) = K(k) は補題 2 から従います。 2n+1(E(an+1, bn+1)− a2n+1T (an+1, bn+1))− 2n(E(an, bn)− a2nT (an, bn)) = 2n−1c2nK(k) を n = 0, . . . , N − 1 に対して加えた後、両辺を −1 倍すると、 −2N(E(a N, bN)− a2NT (aN, bN)) + E(k)− K(k) = − N−1 X n=0 2n−1c2nK(k) (12) となります。ここで、補題 2 と補題 5 を用いました。また、 −2N(E(a N, bN)− a2NT (aN, bN)) = 2N Z π/2 0 a2N p a2 Ncos2θ + b2Nsin 2θ dθ− Z π/2 0 q a2 Ncos2θ + b2Nsin 2θ dθ ! = 2N Z π/2 0 c2Nsin2θ p a2 Ncos2θ + b2Nsin 2θ dθ (13) となります。(13) より、 0≤ −2N(E(aN, bN)− a2NT (aN, bN)) となることが分かります。また、 2N Z π/2 0 c2Nsin2θ p a2 Ncos2θ + b2Nsin 2θ dθ ≤ 2 Nc2 NT (aN, bN) となります。補題 9 と (11) より、 2Nc2NT (aN, bN)≤ c20 2NT (aN, bN) = c20 2NK(k) となります。よって、 −2N(E(a N, bN)− a2NT (aN, bN))≤ c20 2NK(k) となることが分かります。以上より、 0≤ −2N(E(aN, bN)− a2NT (aN, bN))≤ c2 0 2NK(k) となります。この式において、N → ∞ とすると、 −2N(E(a N, bN)− a2NT (aN, bN))→ 0 となることが分かります。よって、(12) において、N → ∞ とすると、 E(k)− K(k) = − ∞ X n=0 2n−1c2nK(k) この式の両辺を K(k) で割ると証明したい式を得ます。 □
定理 3 (ガウスの公式) a = 1, b = 1/√2 としたとき、(9) のように an, bnを定義します。cn = p a2 n− b2nとします。このとき次が成り立ちます。 π = 2AGM 1, 1/ √ 22 1−P∞n=02nc2 n 証明 命題 2 において、k = 1/√2 とおき、両辺に π/2 を掛けると、 π 2 2 E 1/√2 K 1/√2 − 1 ! = π 2K 1/√2 !2 定理 2 と命題 3 より、 π 2 1− ∞ X n=0 2nc2n ! = AGM 1,√1 2 2 よって、定理の式を得ます。 □
