たった一つの整数から円周率を計算する方法
§§０．はじめに．
本稿は以前に書いた「 円周率 π の近似計算 」の要約を兼ねた書き直しである．
あれ以後，私の PC の故障やら何やらもあって思い通りの執筆もままならなかった．
この件を通して私は手塩にかけた愛児を失うが如き悲しみと痛みと残念を経験した．
思えば自らの考案による約１万２千桁に及ぶ一連の円周率の計算は今は亡き愛機への最後の餞（はなむけ）であったと顧みる．
”織り上がる 電子の機（はた）の 円（まろ）き 旋律”（活己）（２０１０年５月１４日）
しかし，お陰で証明や説明の合理性も含めて充分な内省の機会が与えられたことを今ではむしろ感謝している．
あの当時は半信半疑であった幾つかの部分も今ではほぼ完璧に解明されている．
内容は簡潔と平明を心がけ，私にとって自明な事がらの証明や説明は思い切ってこれを省略することとした．
では始めよう．
（ と，まあこんな書き出しで始めてみたのですが，後の方でいろいろ追加やら予定変更がありました．
当初は多角形法的でしかも平方根を計算しなくて済む方法の良いところを強調するつもりでいましたが，後から
逆正接関数の級数展開による近似法などについて自分なりの実行結果や公式の発見などもあって予定を大幅に変更いたしました．
これらの追加や予定変更の部分は私がこの記事を書くときには殆ど私自身の知らなかったことであります．
この記事を書くことによって私自身が良い勉強になり書く前よりはほんの少し円周率の計算に関して賢くなれたことを
告白しておきたいと思います．
インターネットの記事というものは誰かにお読み頂くために書くと言うよりは自分自身のためにこそ
書いてみるべきものであるのかも知れません．(2010'09'16/Thu) ）
§§１ 単位複素弧度回転演算子法による円周率の計算
§１.１ たった一つの任意の整数から始めよう
任意の正の整数を H とし, H に応じて弧度 θ を取り
tan(θ/2) = 1/H (1-1)
であるとする．
このたった一つの定義から以下の式たちが次々に成立する．
sin(θ/2) = 1/√(H2+1) (1-2)
cos(θ/2) = H/√(H2+1) (1-3)
sin(θ) = 2H/(H2+1) (1-4)
cos(θ) = (H2-1)/(H2+1) (1-5)
tan(θ) = 2H/(H2-1) (1-6)
この H を半角分母，θ を単位弧度と呼んでおく．
与えられた単位弧度 θ に対して，商を M, 剰余を δ とすれば
Mθ + δ = π/2 (2)
と書けるが，このような M,δ は一組しかも唯一組しかない．
以下の議論に備え，ここで任意の順序数を m ( m = 1,2,3,… ) として
S(m) = sin(mθ) (3-1)
C(m) = cos(mθ) (3-2)
と書いて良いとする．（これらは全て有理数である．また S(1),C(1) は S,C と書いても良いとしておく．）
m を乗法指数と呼んでおこう．
さてここで，
ζ = cos(θ) + isin(θ) = C + iS (4-1)
なる複素数 ζ を考えよう．（ ζ を単位複素弧度と呼んでおこう．）
しかるに M を適当な順序数, 0 ≦ Mθ ＜ π, 0 ≦ δ ＜ θ
を前提として，ド・モアブルの定理をも考慮に入れれば，
C(m) = Re( C + iS )m ＞ 0 ( m = 1,2,…,M ) (4-2)
C(m) = Re( C + iS )m ＜ 0 ( m = M+1,M+2,…,2M ) (4-3)
Mθ + δ = π/2 ⇔ C(M)C(M+1) ＜ 0 (5-1)
となるような M と δ が一組しかも唯一組だけ常に存在する．
すなわち単位複素弧度 ζ を複素数の乗法二項演算の原理により順次乗じて行けば
有限回で m = M を越える前後で ζm の実部の符号が逆転する．
ゆえに複素数の四則演算の有限回の操作で真の商 M を決定し得る．
少し補足すれば，x を任意の実数として
（１）x ＞ 0 ⇔ sgnum(x) = 1,
（２）x = 0 ⇔ sgnum(x) = 0,
（３）x ＜ 0 ⇔ sgnum(x) = -1
であるような符号関数 sgnum(x) を定義しておけば，
sgnum(cos(mθ)) = sgnum(γ） (γ = π/2 - mθ, m=1,2,3,…M,…2M ) (1-5-2)
が成り立つ．
さらに剰余 δ についても，
sin(δ) ＜ δ ＜ tan(δ), (6-1)
sin(δ) = cos(Mθ) = C(M), (6-2)
tan(δ) = cos(Mθ)/sin(Mθ) ( = cot(Mθ) ) = C(M)/S(M) (6-3)
∴ C(M) ＜ δ ＜ C(M)/S(M) (6-4)
となって剰余 δ も近似できる．
これらに基づいて次の如き二つの不等式が得られ，それらから課題の近似円周率が求められる．
（１）単精度近似不等式
MS ＜ π ＜ (M+1)S/C (7-1)
（２）倍精度近似不等式
MS + C(M) ＜ π ＜ MS/C + C(M)/S(M) (7-2)
これらの不等式を特に終結不等式と呼ぼう．
（ これらの精度についての考察は後に述べよう．）
種々の理由から
ζ = cos(θ) + isin(θ) = C + iS (8-1)
は万能回転演算子とも呼ぶべきものである．
また，乗法指数 m を下げたいときは
ζ-1 = ( cos(θ) + isin(θ) )-1 = C - iS (8-2)
を乗じれば良い．
ζ を前進演算子，ζ-1 を後退演算子と呼ぼう．
ここで述べた円周率の計算方法を単位複素弧度回転演算子法と呼んでおこう．
( ここでは単位の意味を "一つ" ではなく "一つあたり" を意味するものとする．)
以下ではもう少し詳細に細部の関係を調べてみよう．
§１.２ π と H から M を予測する不等式
この節では π と任意の半角分母 H から，これに見合う商 M を求める不等式を導いてみよう．
まず以前に示した商 M と剰余 δ の式を変形して
M + δ/θ = (π/2)/θ (1-1)
を得る．さらに，
sin(θ/2) ＜ θ/2 ＜ tan(θ/2) (2-1)
sin(θ/2) = 1/√(H2+1) (2-2)
tan(θ/2) = 1/H (2-3)
から，
1/√(H2+1) ＜ θ/2 ＜ 1/H (1-2)
を得る．
∴ (π/4)/tan(θ/2) ＜ M + δ/θ ＜ (π/4)/sin(θ/2) (3-1)
∴ (π/4)H ＜ M + δ/θ ＜ (π/4)√(H2+1) (√(H2+1)≒H(1+1/(2H2)) (3-2)
∴ πH/4 ＜ M + δ/θ ＜ πH/4 + π/(8H) (3-3)
∴ 上限 - 下限 = π/(8H) (3-4)
∴ ( 0 ≪ δ/θ ＜ 1, δ/θ ＞ (π/8H) ) ⇒ [πH/4] = M (4-1)
（ [] で直近下位の整数を表す．）
あるいは
( 0 ≒ δ/θ, δ/θ ＜ π/8H ) ⇒ [πH/4] = M-1 (4-2)
となる場合も極めて稀だけれども起こり得る．
これらの式からさらに
既に近似した円周率の下限値 Pa と上限値 Pb が与えられていて
Pa ＜ π ＜ Pb (5-1)
[-log10(Pb-Pa )] ＞ [log10H] (5-2)
であったとすると，次のような近似式が得られる．
PaH/4 ＜ M + δ/θ ＜ (Pb/4)H√(1+1/H2) ( ≒(Pb/4)(H+1/(2H)) ) (6-1)
ゆえに
[PaH/4] = M (6-2-1)
または
[PaH/4] = M-1 (6-2-2)
あるいはもっと簡単に
M ≒ [(Pb+Pa)H/8] (6-2-3)
としてもほぼ 1 の位まで M を予測し得ると考えられる．
（ ここでは単に M を予測するだけなので，この程度の考察で良いとしよう．）
さらに便宜的に半角分母 Hd を
[-log10(Pb-Pa)/2] ＞ log10Hd (7-1)
によって予め定めておくとする．
このとき商の予想値 Mpr を
[(Pb+Pa)Hd/8] = Mpr (7-2)
から予想したとする．
この Hd から定められた単位弧度 θ によって取られた実部 C(m) の符号反転
C(Mtr) ＞ 0 ＞ C(Mtr+1) (7-3)
から真の商 Mtr を決定し，しかも Hd の桁数を単位として単精度と倍精度で円周率が計算できる.
（ この方法を母先商後と呼んでおこう．）
このとき，
1/Hd ≒ θ/2 (7-4-1)
πHd/4 ＞ Mpr (7-4-2)
Mtr + δ/θ = π/(2θ) ( 0＜δ＜θ ) (7-4-3)
から
Mtr + δ/θ ＞ Mpr (7-4-4)
となる．
ここで実数の商 mtr を
mtr = Mtr + δ/θ (7-5-1)
とすれば，
mtr ＞ Mpr (7-5-2)
となる．
計算の実行において，必ずと言って良いほど経験的に,
C(Mpr) ＞ 0 ＞ C(Mpr+1) (7-6-1)
∴ Mtr = Mpr (7-6-2)
であることが確かめられる．
§１.３ π と M から H を予測する不等式
今度は π と任意の M とからこれに見合う H を求める不等式を導いてみよう．
まず
Mθ = π/2-δ (1-1)
sin(θ) ＜ θ ＜ tan(θ) (2-1)
sin(θ) = 2H/(H2+1) (2-2)
tan(θ) = 2H/(H2-1) (2-3)
などから
MH/(H2+1) ＜ π/4-δ/2 ＜ MH/(H2-1) (3-1)
∴ (H2-1)/MH ＜ 4/(π-2δ) ＜ (H2+1)/MH (3-2)
∴ 4M/(π-2δ) - 1/H ＜ H ＜ 4M/(π-2δ) + 1/H (3-3)
∴ H ≒ 4M/π (3-4)
これらの式からさらに
既に近似した円周率の下限値 Pa と上限値 Pb が与えられていて
Pa ＜ π ＜ Pb (4-1)
[-log10(Pb-Pa)] ＞ [log10M] (4-2)
であったならば
4M/Pb ＜ H ＜ 4M/Pa (4-3)
が成り立つと見てよい．
あるいはもっと簡単に
H ≒ [8M/(Pb+Pa)] (4-4)
としてもほぼ 1 の位まで H が予測できる．
（ ここでは単に H を予測するだけなので，この程度の考察で良いとしよう．）
さらに便宜的に剰余 Mpr を
[-log10(Pb-Pa)/2] ＞ log10Mpr (5-1)
によって予想しておいたとする．
この Mpr からさらに半角分母 Hd を
Hd = [8Mpr/(Pb+Pa)] (5-2)
によって予め定めるものする．
この Hd から定められた単位弧度 θ によって取られた実部 C(m) の符号反転
C(Mtr) ＞ 0 ＞ C(Mtr+1) (5-3)
から真の商 Mtr を決定し，しかも Hd の桁数を単位として単精度と倍精度で円周率が計算できる.
（ この方法を商先母後と呼んでおこう．）
このとき，
1/Hd ≒ θ/2 (5-4-1)
Hd ＜ 4Mpr/π (5-4-2)
Mtr + δ/θ = π/(2θ) ( 0＜δ＜θ ) (5-4-3)
から
Mtr + δ/θ ＜ Mpr (5-4-4)
となる．
ここで実数の商 mtr を
mtr = Mtr + δ/θ (7-5-1)
とすれば，
mtr ＜ Mpr (7-5-2)
となる．
計算の実行において，必ずと言って良いほど経験的に,
C(Mpr-1) ＞ 0 ＞ C(Mpr) (7-6-1)
∴ Mtr = Mpr-1 (7-6-2)
であることが確かめられる．
そこでで新たに，
M' = M+1 ( = Mpr ) (7-6-3)
δ' = θ-δ (7-6-4)
となるような M',δ'を採用すれば，
M'θ - δ' = π/2 (7)
であるとも言える．
私は M' を過ぎ越し商，δ' を過ぎ越し剰余と呼んでいる．
結局，商先母後における商の予想値 Mpr は過ぎ越し商であったのである．
§１.３.２ H = fn であるときの M の予想値
( 1 + √2 )n = fn2 - 2gn2 = ± 1 (1-3-2-0)
を基にして
tan(θ/2) = 1/ H = 1/fn (1-3-2-1)
と取ったとすれば，極めて粗い近似として
θ ≒ 2/H = 2/fn (1-3-2-2)
fn ≒ (√2)gn (1-3-2-3)
が成り立つと考えよう．さらに，
Mθ ≒ π/2 (1-3-2-4)
と考えれば
2M/fn = (√2)M/gn ≒ π/2 (1-3-2-5)
∴ M ≒ π/(2(√2))gn (1-3-2-6)
ここで π/(2(√2)) は
正弦波の実効値/最大値 = 1/(√2) (1-3-2-7)
正弦波の平均値/最大値 = 2/π (1-3-2-8)
の比として
正弦波の波形率 = π/(2√2) (1-3-2-9)
と云われ、交流電気回路のテスターなどで比較的安価な平均値型で測定して値だけを実効値で読めるようにするというような応用がある．
しかしながら，
M /gn ≒ 正弦波の波形率 (1-3-2-10)
が成り立つことは意外な発見ではないかと思うし，
さらにもう少し具体的な説明（例えば図形的な説明とか）が無いものかと筆者は常々思案している．
兎に角，このようにして fn, gn を基にして π を近似することにもそれなりの良いところがあるのである．
（この節だけ２０１７年２月４日に加筆）
§１.４ 常用対数を用いた主な数たちの最大桁指数の計算と不等式の精度の計算
一般に任意の正の実数を x とし, 適当な 0 か 正の整数を t として
（１）1 ＜ x ならば
[log10x] = t （ []内で小数点以下を切り捨てる．） (1-1)
によって
（２）0 ＜ x ＜ 1 ならば
-[-log10x] = -t （ []内で小数点以下を切り上げる．） (1-2)
によって x の最大桁指数を t または -t であると定義する．
因みにこの定義によれば，従来の 1 の位の数は最大桁指数が 0 となる．
また 1024 の最大桁指数は 3 となる．
一般に従来の n 桁の整数の最大桁指数は n-1 となる．
さらに 1 未満の小数の場合，例えば 0.00012345 の最大桁指数は -4 となる．
（ 0.00012345 = 1.2345・10-4
∴ [log100.00012345] = [log10(1.2345・10-4)] = [log101.2345+log1010-4]
= [log101.2345-4] = -[4-log101.2345] = -4
このような計算結果にするには
0 ＜ x ＜ 1 では -[-log10x] = -t （ []内で小数点以下を切り上げる．）
とする必要がある．
∵ 3 ＜ 4-log101.2345 ＜ 4 ）
それではこれを用いて前述の主な数たちの最大桁指数を求めてみよう．
まず，以後の議論では常に
[log10H] = t (2)
を共通の前提とする．
すると前節で既に示したように
（１‐１）M ≒ (π/4)H (3-1-1)
から
[log10M] = [log10(π/4)H] = [log10(π/4)+log10H] = t (3-1-2)
（１‐２）tan(θ/2) = 1/H (3-2-1)
から
[log10(tan(θ/2))] = [log101-log10H] = -[log10H] = -t (3-2-2)
（１‐３）sin(θ/2) = 1/√(H2+1) (3-3-1)
から
[log10(sin(θ/2))] = [log10(1/√(H2+1))] = [log10(1/H√(1+1/H2))]
= -[log10H√(1+1/H2)] = -[log10H+log10√(1+1/H2)]
= -[log10H] = -t (3-3-2)
（１‐４）cos(θ/2) = H/√(H2+1) (3-4-1)
から
[log10(cos(θ/2))] = [log10(H/√(H2+1))] = [log10(1/√(1+1/H2)]
= [-log10√(1+1/H2)] = -[(1/2)log10(1+1/H2)] ≒ 0 (3-4-2)
（ 実際には cos(θ/2) = 0.999999… のようになっている．ゆえに最大桁指数はほぼ 0 となる．（微妙））
（１‐５）tan(θ) = 2H/(H2-1) (3-5-1)
から
[log10(tan(θ))] = [log10(2H/(H2-1)] = [log10(2/(H-1/H)]
= [log102-log10(H(1-1/H2))] = [log102-log10H-log10(1-1/H2)]
= -[log10H+log10(1-1/H2)-log102] = -[log10H] = -t (3-5-2)
（１‐６）sin(θ) = 2H/(H2+1) (3-6-1)
から
[log10(sin(θ))] = [log10(2H/(H2+1))] = [log10(2/(H+1/H))] = [log102-log10H(1+1/H2)]
= [log102-log10H-log10(1+1/H2)] = -[log10H+log10(1+1/H2)-log102] = -[log10H] = -t (3-6-2)
（１‐７‐１）cos(θ) = (H2-1)/(H2+1) (3-7-1)
から
[log10(cos(θ))] = [log10(H2-1)/(H2+1)]
= [log10H2(1-1/H2)-log10H2(1+1/H2)]
= [2log10H+log10(1-1/H2)-2log10H-log10(1+1/H2)] ≒ 0 (3-7-2)
（１－７－２）因みに cos(θ) で 9 が始めから何桁並んでいるかも調べておこう．
[log10(1-cos(θ))] = [log10(1-(H2-1)/(H2+1))]
= [log10(2/(H2+1))] = [log102-log10(H2+1)] = [log102-log10H2(1+1/H2)]
= [log102-2log10H-log10(1+1/H2)] = -[2log10H+log10(1+1/H2)-log102]
= -2[log10H] = -2t (3-7-3)
ゆえに cos(θ) の始めの 2t 桁ほどは 9 ばかりが並んでいる．
次にそれぞれの不等式についてもその近似の誤差の最大桁指数や精度を定義しよう．
一般に x を任意の実数, A, B を適当な有理数たちとして近似不等式
A ＜ x ＜ B, ( 0 ＜ B - A ＜ 1 ) (4-1)
が与えられているとする．
このとき 0 を含む適当な正の整数 t が存在して
[-log10(B-A)] = t (4-2)
となるはずである．
このとき 4-1 の近似不等式の誤差の最大桁指数は -t であると言うべきである．
さらに 4-1 の近似の精度は小数点以下 t-1 桁であると言うべきである．
なぜなら 4-1 において A と B は小数点以下 -t+1 の桁指数までは一致しており，
-t の桁指数で不一致になっているはずだからである．
この事実を簡単に 4-1 は t-1 桁の精度を持つと言って良いとしよう．
あるいは x は（小数点以下）t-1 桁目まで求まっていると言っても良いと考えられる．
それではこれらの定義に従って，先述の代表的な不等式たちの精度を計算してみよう．
ここでも先ず
[log10H] = t (5-1)
を前提とする．
（２‐１）sin(θ/2) ＜ θ/2 ＜ tan(θ/2) (5-2-1)
sin(θ/2) = 1/√(H2+1) ≒ 1/(H+1/(2H)) (5-2-2)
tan(θ/2) = 1/H (5-2-3)
から
2H/(2H2+1) ＜ θ/2 ＜ 1/H (5-2-4)
このとき
1/H - 2H/(2H2+1) = (2H2+1-2H2)/(H(2H2+1))
= 1/(2H3+H) = 1/(H3(2+1/H2)) (5-2-5)
∴ [-log10(1/(H3(2+1/H2)))] = [log10(H3(2+1/H2))]
= [log10H3+log10(2+1/H2)] = [3log10H +log10(2+1/H2)] = 3t (5-2-6)
となるから 5-2-4 で近似した場合に θ/2 は（小数点以下）3t-1 桁の精度で求まっている．
ここから θ/2 は（ t を単位として ）三倍精度で求まっていると言えよう．
（ 結局のところ，
0 ＜ δ ＜ θ ⇒ tan(δ) - sin(δ) ＜ tan(θ) - sin(θ)
となるので，むしろ剰余 δ の方が θ より正確に求められるというべきである.
さらに
sin(δ) = cos(Mθ), tan(δ) = cot(Mθ) = cos(Mθ)/sin(Mθ)
であることに注意すれば，
cos(θ), sin(θ) から cos(Mθ), sin(Mθ) を求めたならばアルキメデスの方法のように
多重平方根を開く回りくどさと比べて計算がうんと簡単で速くなると言うべきである．
この事実はアルキメデスのように正多角形を用いる方法が準正多角形を用いるここでの方法と比べてより正確であろう
というような迷信を払拭する事実である．（発見！！）
そして多角形を用いた円周率の近似計算法の中で明確には指摘されて来なかった事実ではないのかと私は言いたいのである．
すでにアルキメデスと同時代のピタゴラスが任意の奇数 H から三つの整数
(H2-1)/2, H, (H2+1)/2
の間で三平方の定理が成り立つことを発見しているのであるからアルキメデスほどの実力があるならば
私が述べたような方法も見破れたかも知れないのである．
もし，そんなことにでも成っていたら多角形で近似するような円周率の計算の歴史も，
もう少し実用的なものに変わっていたかも知れないのではないかと私は思うのである．
しかし，その為には少なくとも負の数の概念のようなものがなければ，
cos(Mθ) ＞ 0, cos((M+1)θ) ＜ 0
から M, cos(Mθ), sin(Mθ) の値を決定することは困難であったろうとも思われる．
要約的に，これらの方法を代数に例えるならば，
アルキメデスの方法は二次分解を繰り返して行くことに相当するので平方根が必要となる．
一方，私のやり方は式を展開して行くことに相当するので冪根は全く不要で有理計算だけで足りるのである．）
同じようなことだが
（２‐２）sin(θ) ＜ θ ＜ tan(θ) (5-3-1)
sin(θ) = 2H/(H2+1) (5-3-2)
tan(θ) = 2H/(H2-1) (5-3-3)
から
2H/(H2+1) ＜ θ ＜ 2H/(H2-1) (5-3-4)
このとき
2H/(H2-1) - 2H/(H2+1) = 4H/(H4-1) (5-3-5)
∴ [-log10(4H/(H4-1))]
= [-log10(4H) + log10(H4-1)] = 3[log10H] = 3t (5-3-6)
となるから 5-3-4 で近似した場合に θ は（小数点以下）3t-1 桁の精度で求まっている．
ここから θ も θ/2 と同様に三倍精度で求まっていると言えよう．
（２‐３）Msin(θ) ＜ π/2 ＜ (M+1)tan(θ) (5-4-1)
M ≒ (π/4)H (5-4-2)
から
(M+1)tan(θ) - Msin(θ) ≒ (π/4)H( tan(θ) - sin(θ) ) + tan(θ) (5-4-3)
ところが（１‐２）の結果から
[-log10(π/4)H(tan(θ) - sin(θ))] = [-log10H + 3log10H]
= 2[log10H] = 2t (5-4-4)
一方
[-log10tan(θ)] = t (5-4-5)
ゆえに 5-4-1 による π の近似の精度は t-1 桁であって，単精度となる．
（２‐４）Msin(θ) + cos(Mθ) ＜ π/2 ＜ Mtan(θ) + cot(Mθ) (5-5-1)
M ≒ (π/4)H (5-5-2)
0 ＜ cot(Mθ) - cos(Mθ) ＜ tan(θ) - sin(θ) (5-5-3)
から
[-log10((Mtan(θ)+cot(Mθ)) - (Msin(θ)+cos(Mθ)))]
≒ [-log10(M+1)(tan(θ)-sin(θ))] = [-log10(M+1) - log10(tan(θ)-sin(θ))]
= [-log10H + 3log10H] = 2t (5-5-4)
ゆえに 5-5-1 による π の近似の精度は 2t-1 桁であって，倍精度となる．
（ これ以降の議論では特に正確に 1 桁の違いに拘る必要がない場合も多い．
このような場合，上限と下限の差の値の最大桁指数が -t であれば精度も同様に t (桁)であると言っても良いとする．）
§１.５ 倍角公式と加法定理の組み合わせによる演算回数の圧縮
ここで例えば，ある半角分母を H 商を M として
M = 21 (1-1)
であったと仮定しよう．
すると，既に説明したように，
C = (H2-1)/(H2+1) (2-1)
S = 2H/(H2+1) (2-2)
ζ = C + iS (2-3)
とすれば，
Re(ζm) ＞ 0 ( m = 1,2,3,…,21) (3-1)
Re(ζ22) ＜ 0 (3-2)
となるはずである．
さて，ここで倍角公式
cos(2mθ) + isin(2mθ) = (cos(mθ)+isin(mθ))2 (4-1)
を使って演算回数を軽減させることを考える．
そこで複素数の乗法二項演算が何回起こるかを実際に追ってみよう．
M = 21 = 16 + 4 + 1 = 1×24 + 1×22 + 1×20 (4-2)
から
（１回目）(cos(θ)+isin(θ))2 = cos(21θ)+isin(21θ) (5-1)
で倍角公式を１回実行．
（２回目）(cos(21θ)+isin(21θ))2 = cos(22θ)+isin(22θ) (5-2)
で倍角公式を１回実行．
（３回目）(cos(22θ)+isin(22θ))(cos(θ)+isin(θ)) = (cos(5θ)+isin(5θ) ) (5-3)
で加法定理を１回実行．
（４回目）(cos(22θ)+isin(22θ))2 = cos(23θ)+isin(23θ) (5-4)
で倍角公式を１回実行．
（５回目）(cos(23θ)+isin(23θ))2 = cos(24θ)+isin(24θ) (5-5)
で倍角公式を１回実行．
（６回目）(cos(24θ)+isin(24θ))(cos(5θ)+isin(5θ)) = cos(21θ)+isin(21θ) (5-6)
で加法定理を１回実行．
（７回目）cos(21θ)+isin(21θ))(cos(θ)+isin(θ)) = cos(22θ)+isin(22θ) (5-7)
で加法定理を１回実行．このとき実部の符号が始めて負となることによって終結する．
このような倍角公式と加法定理の組み合わせでは本来 22 回の複素数の乗法二項演算が僅か７回で済んでいる．
これをもう少し一般化して述べておこう．
M を二進数で表したとき，0 でない最大桁指数が t であれば
MaxBit(M)2 = t (6-1)
で表して良いとする．
今の例では，
M = 21 = 16 + 4 + 1 = 1×24 + 1×22 + 1×20 (4-2)
なので
MaxBit(M)2 = MaxBit(21)2 = 4 (6-2)
となる．
そして倍角公式の実行回数を Nt とすれば
Nt = MaxBit(M)2 = 4 (6-3)
さらに M を二進数で表したとき，1 である桁の個数が n であれば
NumBit(M)2 = n (7-1)
で表して良いとする．
今の例では，
NumBit(M)2 = NumBit(21)2 = 3 (7-2)
となる．
そして加法定理の実行回数を Na とし，最後の実部の符号反転に達する分も勘定に入れれば
Na = NumBit(M)2 = 3 (7-3)
それで倍角公式と加法定理を組み合わせた場合の複素数の乗法二項演算の実行回数を Nta とすれば，
Nta = Nt + Na = MaxBit(M)2 + NumBit(M)2 (8-1)
今の例では，
Nta = MaxBit(21)2 + NumBit(21)2 = 4 + 3 = 7 (8-2)
となる．
このとき，
η = Nta/M (9-1)
を倍角圧縮比と呼んでおこう．
今の例では，
η = 7/21 (9-2)
となる．（表示は分数のままで良いとする．）
ところが，一般に，
1 ≦ NumBit(M)2 ≦ MaxBit(M)2+1 (9-3)
また M = 1024 = 210 ならば
MaxBit(1024) = 10 (9-4)
そこで
M ≒ 1024k (9-5)
ならば，
MaxBit(M) = MaxBit(1024k) = MaxBit(210k) = 10k (9-6)
このとき平均的に
NumBit(M)2 ≒ MaxBit(M)/2 = 5k (9-7)
であるとし，圧縮比を
ηmean ≒ (1+1/2)MaxBit(M)/M = 15k/1024k (9-8-1)
と考え，平均倍角圧縮比と呼んでおこう．
同様に，
ηmin ≒ (1+1)MaxBit(M)/M = 20k/1024k (9-8-2)
を最小倍角圧縮比，
ηmax ≒ MaxBit(M)/M = 10k/1024k (9-8-3)
を最大倍角圧縮比と呼んでおこう．
これは複素数どうしの二項演算の回数において 1024k 回の二項演算が僅か 10k 回の二項演算に縮むということである．
結局のところ, 商 M が丁度 2 の冪であるとき η が最大倍角圧縮比となると考えて良い．
特に M = 2k であったとき
ηmax = k/2k (9-8-4)
となる．
このような商 M を二冪商と呼ぼう．
M が二冪商に取られているとき複素数の乗法二項演算の回数はこの付近の商の内では最小となる．
さらに一般的に考えて見よう．
x, y を適当な実数として，
M = 10x = 2y (10-1)
であったとする．これらの常用対数を取れば，
log10M = log1010x = log102y
∴ log10M = x = ylog102 ( log102 ≒ 0.30103 ) (10-2)
ここでも t を正の整数として，
[log10M] = t (10-3)
であったとすれば，倍角公式と加法定理の組み合わせによる複素数の乗法二項演算の回数 Nta は
Nta ≒ y = log10M/log102 = t/log102 (10-4)
となる．
∴ η ≒ Nta/M = log10M/(Mlog102) = t/(Mlog102) (10-5)
このηを理論圧縮比と呼んでおこう．
一例を述べれば約 12000 桁の円周率の計算では M は約 6000 桁の整数となるが，これを，
M ≒ 10x ≒ 2y (10-6-1)
x ≒ 6000 (10-6-2)
と見なせば，
M ≒ 106000 = 10002000 ≒ 10242000 = 220000
∴ Nta ≒ y = 20000 (10-6-3)
となる．
これは複素数の乗法二項演算の回数の意味において，
6000 桁の整数で表された回数の二項演算が僅か 5 桁の整数で表された回数の二項演算に縮むということである．
一方，a,b,c を実数として，
( cos(a) + isin(a) )( cos(b) + isin(b) ) = cos(c) + isin(c) (11-1)
cos(c) = cos(a)cos(b) - sin(a)sin(b) (11-2)
sin(c) = sin(a)cos(b) + cos(a)sin(b) (11-3)
c = a + b (11-4)
の複素数の乗法二項演算において１０進１桁どうしの乗法二項演算の回数を Nma とし，これを求めて見よう．
上記複素数の実部と虚部の計算領域が t に一次比例する桁数であったとすれば, 比例定数を k として，
Nma ≒ kt2 (11-5)
と考えて良い．
それで１回の近似計算で終結不等式に至る一連の演算における１０進１桁どうしの乗法二項演算の回数を Nall とすれば，
Nall ≒ NtaNma = kt3 (11-6)
ゆえに
t' = 2t (12-1)
ならば
Nall'/Nall ≒ (t'/t)3 = 8 (12-2)
すなわち分母 H または商 M の桁数が二倍になれば乗法二項演算の回数は約８倍となる．
その結果，計算の所要時間も約８倍となる．
私の最近の実行結果を例に取れば，この記事を書いている今の PC
（ 今ではもう故障してしまった DELL 社製デスクトップ PC の後継機として現プロバイダ経由で入手した HP 社製ノート PC ）
で Java-BigDecimal に依れば M = 2N として
N = 2560 を入力実引数とした円周率 1541 桁の近似計算では約 30 秒，
N = 5120 を入力実引数とした円周率 3083 桁の近似計算では約 3 分 48 秒，
N = 10240 を入力実引数とした円周率 6165 桁の近似計算では約 29 分 2 秒，
N = 20480 を入力実引数とした円周率 12330 桁の近似計算では約 3 時間 58 分 9 秒，
となって精度が約 2 倍になると処理時間が約 8 倍となっている．
更に次の極めて正確な計算が成り立っている．
2log1022560 = 2・2560・log102 = 1541.273… ≒ 1541
2log1025120 = 2・5120・log102 = 3082.547… ≒ 3083
2log10210240 = 2・10240・log102 = 6165.094… ≒ 6165
2log10220480 = 2・20480・log102 = 12330.188… ≒ 12330
すなわち M = 2N, 近似される円周率の精度を前回を Told 今回を Tnew として
N ≒ Told/log102 (12-3-1)
Tnew ≒ 2log10M = 2log102N = 2Nlog102 (12-3-2)
Tnew ≒ 2Told (12-3-3)
が極めて正確に成り立っている．（ これは近似の検算に使える．）
（ 以下をクリックすると最近の実行結果のテキストファイルが開けます．データが横に広がって見にくいので
(InternetExplorerのツールバー) ⇒ [表示] ⇒ [ソース(C)] ⇒ (別窓でファイルを表示) を使うと良いでしょう．）
N=2560,T=1541
N=5120,T=3083
N=10240,T=6165
N=20480,T=12330
（ ※ このプログラムでは標準出力をコマンドプロンプト画面ではなくテキストファイルにリダイレクトしたので
実行時間が前者と比較して約 1/3 程度早くなっている．）
これら一連の計算における方針は
（１）前回の近似円周率の近似精度を Told, 冪指数を `N とし, 今回の冪指数を N, 商を M とし
N = 2`N ＜≒ Told/log102 (13-1-1)
（ 些細なことだが自然対数しか使えないプログラム言語では log102 = logε2/logε10 を用いよ．ε は自然対数の底 ）
M = 2N (13-1-2)
となるように取る．
（２）前回の円周率の近似下限値を Pa, 上限値を Pb として半角分母 H を
H = [8M/(Pa+Pb)] ( Pa＜ π ＜Pb ) (13-2)
となるように取る．
（３）単位複素弧度 ζ を
ζ = C + iS = (H2-1)/(H2+1) + i2H/(H2+1) (13-3-1)
となるように取る．
（４）ζ2m = C(2m) + iS(2m) = ( C(m) + iS(m))2 (13-4)
を 2m = 21,22,23,…,2N まで合計 N 回実行する．
このとき倍角公式から
C(2m) = C(m)2 - S(m)2 (13-5-1)
S(2m) = 2C(m)S(m) (13-5-2)
となるように取れば良い．
ここまでの最後の商は m = 2N となる．
（５－１）C(m) が正であったならば ζ を 1 回乗じて乗法指数 m を １ だけ加算し C(m) が負になるまで繰り返す．
（５－２）C(m) が負であったならば ζ-1 を 1 回乗じて乗法指数 m を １ だけ減算し C(m) が正になるまで繰り返す．
ここで
ζ-1 = C - iS (13-3-2)
として良い．
（６）C(m) の符号が正である m の最大値を真の商 M であるとし M を決定する．
MS + C(M) ＜ π/2 ＜ MS/C + C(M)/S(M) (13-6-1)
から新しい近似下限値 Pa', 上限値 Pb' を得る．
このとき上式の近似円周率を Pinew, 誤差を Er, 精度を Tnew とすれば，
Pa' = 2( MS + C(M) ) (13-6-2)
Pb' = 2( MS/C + C(M)/S(M) ) (13-6-3)
Pinew = (Pb'+Pa')/2 (13-6-4)
Er = (Pb'-Pa')/2 (13-6-5)
Tnew = [log10(Pb'-Pa')/2] ≒ 2Told (13-6-6)
となるはずである．（ 以後 Told を単に t と言おう．）
（７）私は C や S の少数点以下の計算領域（スケール）を 3t+3 桁（以上）に取った．
こうすると t 桁の整数の商との掛け算で末尾の 3 桁を丸め用に使っても少なくとも 2t 桁の精度が保証されるからである．
C(2),S(2) 以降は 2t+3 桁に取っても良さそうに思えるが経験的には不正な結果となる．
これらも 3t+3 桁（以上）に取るのが安全である．
（ この理由を考えるに，例えば S(2) ではまだ少数点以下 t 桁目まではほぼ 0 ばかりが並んでいるので，
スケールを 2t+3 桁に取ると先頭の 0 たちを除いた後の部分が t+3 桁ほどしかなく，
正確な計算に必要な 2t+3 桁（以上）を満足しないからである．
必要な精度を満たすためには 0 でない最大の桁から 2t+3 桁（以上）のスケールに取るべきである．
因みに S は始めの t 桁ほどは 0 ばかりが並んでいる．C は始めの 2t 桁ほどは 9 ばかりが並んでいる．
なおスケールを必要以上に大きく取りすぎると処理時間が長くなる，しかし計算結果の正しさには影響しない．）
いずれにせよ毎回の不等式が独立しているので次回の倍精度の計算結果からか，
または別の独立した同精度以上の計算結果を照合することで計算間違いの有無を確認できる．
この一連の円周率近似法を”単位複素弧度回転演算子による二冪商二倍精度追跡法”と呼んでおこう．
さらに一般に商が二冪商でない場合であっても倍角公式と加法定理によって演算回数を減らしている場合を
”単位複素弧度回転演算子による倍角加法定理圧縮法”と呼んでおこう．
以後では後者の意味を表す用語として単に”単位複素弧度回転演算子法”と呼んでおくことにしよう．
この方法は
理解の容易,平方根近似の不要,微積分の不要,無限概念の不採用,有限唯一性原理の採用,不等号検査の同時,毎回の不等式の独立,
検算の自己完結,プログラムの適正な動作確認の容易,計算能率の高速
等の理由から特殊な課題の習得を目的としない簡易な円周率の近似においては他の方法のどれよりもより実用的であると思える．
そのようなわけで私はこの方法を世界で一番良い円周率の計算方法であると自負している．
この自信は初めてこれと類似の方法を思いついた頃から日増しに強くなりこそすれ少しも揺らいでいない．
特に，この方法と同義な原理はアルキメデスの方法や，
オイラーやマチンの方法といったような他の円周率の計算方法の証明でさえより簡明なものと成し得る．
加えてこの方法は初等幾何学の僅かな前提だけから導き得るので誰もが心得ておくべき基本原理的な意味合いが強いと言えよう．
次章でそのことをもう少し詳しく述べよう．
§§２ 単位複素弧度回転演算子法の優越性
§２.１ 各種の円周率計算法における二項演算的同型性
任意の実数を θ とし，次のような写像 φ によって複素数 ζ を対応させる．
ζ = φ(θ) = εiθ ( εは自然対数の底 ) (1-1-1)
上式はまたオイラーの公式により
εiθ = cos(θ) + isin(θ) (1-1-2)
であるとも言える．
ゆえに，このような ζ には以前の単位複素弧度回転演算子も含まれる．
この写像 φ は実数 θ の加法群から複素数 ζ の乗法群への群同型な写像であると見なす事ができる．
ところで，一見異なって見える様々な円周率の近似計算法も大きくは次の二種類に分けることができる．
すなわち
（１）単位半径の複素円周上の多角形の頂点による近似法
これにはアルキメデスの方法，関孝和の方法，単位複素弧度回転演算子法などが含まれる．
（ 以後これらを総称して多角形法と呼ぼう．）
（２）弧度の無限級数の有限部分和による近似法
これにはグレゴリーの公式，シャープの公式，オイラーの公式，マチンの公式などが含まれる．
（ 以後これらを総称して級数法と呼ぼう．）
しかるに写像 φ において，指数 θ を（２）に対応させ，複素数 ζ を（１）に対応させるならば，
あるいはもう少し丁寧に述べれば，
φ(θ1+θ2) ⇔ φ(θ1)φ(θ2) (1-1-3)
が成り立つので，左辺の括弧内の加法を級数法に対応させ，右辺の乗法を多角形法に対応させるならば，
これらの大きく分けられた二種類の近似法は二項演算的には同型であると見なし得る．
少なくともこの同型性から
◇”多角形法における頂点と級数法における項とは互いに同義的である．” (2)
ことだけは示されたと言えよう．
（ あるいは関孝和の方法は多角形法であると同時に級数法であるとも考えられる，このことからもこの主張の正当性が伺える．）
このようにして，一先ず凡そ全ての円周率の近似法を同列に並べて置いてから，
理論的と実用的の両面からそれぞれの方法の長所と短所を比較して見るのも一興であろうと思われる．
§２.２ アルキメデスの方法と単位複素弧度回転演算子法との比較
ここでは円周率の計算に関してアルキメデスの方法を応用した場合と単位複素弧度回転演算子法を応用した場合について
数学的な類似点と相違点についてを簡単にまとめておこう．
さらに自作のプログラムによる計算結果を用いて主に性能の比較と要求精度に必要な計算領域（スケール）を考察してみよう．
参考として，ソースプログラム (Javaファイル）と実行結果の出力ファイルも公表しておいた．
使用言語は Java とし，使用数値は主に BigDecimal 型とした．
（１）アルキメデスの方法と単位複素弧度回転演算子法の数学的類似点と相違点
先ず両者は
Mθ + δ = π/2 (1-1-1)
であるとしたとき
Msin(θ) + sin(δ) ＜ π/2 ＜ Mtan(θ) + tan(δ) (1-1-2)
のような不等式を用いて π の値を近似するという点では類似していると言えよう．
ただ異なっている点は
アルキメデスの方法では
δ = 0 (2-1)
M = 3・2n (n=1,2,3…) (1-2-2)
θ = (π/6)/2n (1-2-3)
主に半角公式
cos((π/6)/2m) = √((1+cos((π/6)/2m-1))/2) (1≦m≦n) (1-2-4)
などを繰り返し用いることになるだろうということである．
一方，単位複素弧度回転演算子法では
δ ≠ 0 (1-3-1)
M ≒ (π/4)H (1-3-2)
sin(θ) = 2H/(H2+1) (1-3-3-1)
cos(θ) = (H2-1)/(H2+1) (1-3-3-2)
tan(θ) = 2H/(H2-1) (1-3-3-3)
主に倍角公式
cos(2nθ) + isin(2nθ) = (cos(2n-1θ) + isin(2n-1θ))2 (1-3-4)
などを繰り返し用いることになるだろうということである．
（２）アルキメデスの方法の中で平方根の計算について，二分法とニュートン法の性能の比較
要求精度が 100 桁で平方根を二分法で計算した場合
所要時間 0 分 19 秒 390 ミリ秒
要求精度が 100 桁で平方根をニュートン法で計算した場合
所要時間 0 分 0 秒 156 ミリ秒
要求精度が 200 桁で平方根を二分法で計算した場合
所要時間 5 分 50 秒 844 ミリ秒
要求精度が 200 桁で平方根をニュートン法で計算した場合
所要時間 0 分 0 秒 874 ミリ秒
このようにニュートン法の方が二分法よりも圧倒的に速いので以後はニュートン法を用いるものとする．
（３）単位複素弧度回転演算子法とアルキメデスの方法における実用的な観点からの三つの課題
これらの方法を応用して主に PC 上で実行できるような高級言語（ C や JAVA 等 ）による計算用
アプリケーションプログラムを組み立てる場合において，
少なくとも次のような三つの課題を解決しておくべきであると思われる，
[1] 要求精度が T 桁であるとき，商 M は如何ほどか？
[2] 要求精度が T 桁であるとき，必要な倍角公式または半角公式の回数は如何ほどか？
[3] 要求精度が T 桁であるとき，必要な計算領域（スケール）は如何ほどか？
これらの内のいづれにおいても
単位複素弧度回転演算子法でしかも M が二冪商に取られている場合を初めに考察しておくと
アルキメデスの方法における同種の問題も簡単に解決してしまうのである．
それでは上の三つの問題を順に解いて行こう．
[1] 要求精度が T 桁であるとき，商 M は如何ほどか？
既に単位複素弧度回転演算子法の説明で述べたように
log(M)≒log(H)≒t (3-1-1)
T ≒ 2t (3-1-2)
が成り立つ．
ゆえに先ず単位複素弧度回転演算子法において
log(M)≒T/2 (3-1-3-1)
が成り立つ．
しかも単位複素弧度回転演算子法は任意の商 M においても成り立つので
アルキメデスの方法にも近似的に応用できると考えられる．
ゆえにアルキメデスの方法においても
log(M)≒T/2 (3-1-3-2)
が成り立つと考えて良い，
[2] 要求精度が T 桁であるとき，必要な倍角公式または半角公式の回数は如何ほどか？
単位複素弧度回転演算子法で M を二冪商に取り，必要な倍角公式の回数を N とすれば
M = 2N (3-2-1-1)
となる．
さらに 3-1-3-1 から
N ≒ T/(2log2)≒ T/(2*0.301…)≒5T/3≒1.666…T≒1.7T (3-2-1-2)
ここではアルキメデスの方法の正多角形を四分の一円内で考えることにし
M = 3*2N (3-2-2-1)
としておけば N が必要な半角公式の回数で，必要な平方根の計算回数より一回だけ少ない．
3-1-3-2 から
N ≒ T/(2log2)≒ T/(2*0.301…)≒5T/3≒1.666…T≒1.7T (3-2-2-2)
3-2-1-2, 3-2-2-2 から分かることは要求精度 T が同じである限り
単位複素弧度回転演算子法で M を二冪商に取った場合の倍角公式の回数と
アルキメデスの方法の半角公式の回数はほとんど同じであるということである．
[3] 要求精度が T 桁であるとき，必要な計算領域（スケール）は如何ほどか？
この問題は極めて実用的な問題でありながら今までの二つの問題と比較すると格段の難問であろうと思われる．
しかしここでも単位複素弧度回転演算子法でしかも二冪商の場合が比較的に綺麗で楽に解け
そこからまたアルキメデスの方法についての同種の問題も意外にあっけなく解けてしまうのである．
では，やってみよう．
以前にも触れたように，
log(sin(θ)) = log(2H/(H2+1))
≒ log(H)-2log(H) = -log(H)≒ -t (3-3-1-1)
log(tan(θ)) = log(2H/(H2-1))
≒ log(H)-2log(H) = -log(H)≒ -t (3-3-1-2)
となるが，これら 3-3-1-1, 3-3-1-2 が意味していることは
単位複素弧度回転演算子法における sin(θ) や tan(θ) は
少数点以下の初めの t-1 桁ほどは 0 ばかりが並んでいて
少数点以下 t 桁目にして初めて 0 で無い数が現れるということである．
さらに
log(tan(θ)-sin(θ)) = log(2H/(H2-1)-2H/(H2+1))
= log(4H/(H4-1))≒ log(H) - 4log(H) = -3t (3-3-1-3)
この 3-3-1-3 が意味していることは不等式
sin(θ) ＜ θ ＜ tan(θ) (3-3-1-4)
は少数点以下 3t 桁目で初めて不一致になるということである．
これらのことから単位複素弧度回転演算子法の終結不等式
Msin(θ) + cos(Mθ) ＜ π/2 ＜ Mtan(θ) + cot(Mθ) (3-3-1-5-1)
において
log( Msin(θ) + cos(Mθ) )≒ log(Msin(θ))
= log(M) + log(2H/(H2+1))≒ log(M) + log(H) - 2log(H)
= t + t - 2t = 0 (3-3-1-5-2)
および
log( Mtan(θ) + cot(Mθ) )≒ log(Mtan(θ))
= log(M) + log(2H/(H2-1))≒ log(M) + log(H) - 2log(H)
= t + t - 2t = 0 (3-3-1-5-3)
から終結不等式 3-3-1-5-1 の下限と上限は 1 の位から始まる．
さらに
log( Mtan(θ) + cot(Mθ) - Msin(θ) - cos(Mθ) )
≒log( M(tan(θ) - sin(θ)) ) = log(M) + log( tan(θ) - sin(θ) )
= log(M) + log( 2H/(H2-1) - 2H/(H2+1 )
= log(M) + log(4H/(H4-1))
≒ log(M) + log(H) - 4log(H)
= t + t - 4t = -2t (3-3-1-5-4)
ゆえに終結不等式 3-3-1-5-1 は少数点以下 2t-1 桁まで下限と上限が一致するはずである．
ゆえに求むる円周率もまた 2t-1 桁の精度となるであろう．
ここで要求精度を T として
t ≒ 0.5T (3-3-1-6)
であることを再度注意しておきたい．
結局のところ，単位複素弧度回転演算子法では初めの sin(θ) や tan(θ) の計算領域（スケール）は
少なくとも 1.5T 桁以上に取るべきであると考えられるが，
その内で sin(θ) や tan(θ) の少数点以下の初めの 0.5T 桁ほどは 0 しか入っていないので
実質においては 1.5T 桁の内の前方の 0 たちを除いた残りの約 T 桁ほどだけが
計算に使われているにすぎないとも考えられる．
そして最後の要求精度に直接係わる部分の計算領域では T+3 桁以上もあれば良いと考えられる．
この結果，単位複素弧度回転演算子法では
計算の一番初めの 1.5T 桁が期間全体の最長スケールであると考えられる．
それでスケールを初めから終わりまで 1.5T 桁で固定しておいても T 桁の要求精度が満たされはするが，
すでに [2] で明らかなように M が二冪商に取られているときは
倍角公式は約 1.7T (≒T/(2log2) ) 回ほども実行される．
このとき期間の終わりに達する間に 0.5T 桁縮まれば
1.5T - 0.5T = T (3-4)
となって T 桁の要求精度が満たされるだけでなく処理速度も向上すると考えられる．．
なお期間全体で縮める合計は 0.5T 桁よりも若干少な目でも良いが均等に縮めるのが良いと考えられる．
そこで，倍角公式 4 回毎に 1 桁づつ，または 17 回毎に 5 桁づつスケールを縮めるような
可変スケール法で実行しても良いと考えられる．
私の実行結果では 1.5T 桁の固定スケール法と比較して
二つの内のどちらの可変スケール法でも同じ要求精度を満たして，
しかも約 1/3 ほど処理速度が向上するようである．
（ 二つの内では 4 回で 1 桁よりも 17 回で 5 桁の方がほんの少しより速くなっている．）
言い替えれば二冪商の場合は倍角公式 17 回毎に 5 桁づつ計算領域が余ってくる，
期間全体の合計では約 0.5T 桁も余るのだとも考えられるのである．
そしてこの事実がまたアルキメデスの方法のスケール考察のヒントとなるのである．
（４）アルキメデスのジレンマ
アルキメデスの方法における計算の出発点は常に下限側では
sin(π/6) (= sin30°) = 1/2 (4-1-1)
上限側では
tan(π/6) (= tan30°) = √3/3 (4-1-2)
となるが，
これらは本来は無限に長く続いていた無理数であったが
計算の便宜から s 桁の有限な計算領域に取って計算を実行するのである．
すると s が短かすぎれば要求精度 T 桁を満たさず，
逆に長すぎれば要求精度 T 桁を満たしはしても計算時間が無用に長く掛かりすぎてしまうであろう．
合理的な計算領域（スケール）は如何ほどであり，またその根拠は如何なるものか？
アルキメデスの方法における半角公式は
単位複素弧度回転演算子法で M が二冪商に取られている場合の倍角公式の
逆の演算にあたると考えられる．
するとアルキメデスの方法では半角公式を一回実行するごとに倍角公式とは逆に
約 17 回毎に 5 桁づつスケールを長くしていかなければならないと考えられる．
期間全体では合計で約 0.5T 桁だけ最後に必要な桁数より長くなると考えられる．
しかもアルキメデスの最後に計算された中心角(≒π/6/2N) は
同じ要求精度を T 桁として
単位複素弧度回転演算子法で M が二冪商に取られている場合の初めの中心角(≒2tan-1(1/H))
とほぼ同じ程度の大きさであると考えて良い．
するとアルキメデスの方法の最後の中心角における必要な計算領域は約 1.5T 桁となる．
この結果アルキメデスの方法における合理的な計算領域の候補の一つとして，
これを Sa とすれば
Sa = 1.5T + 0.5T + α = 2T + α (4-2-1)
が挙げられる．
この Sa を経験的最小スケールと呼んでおこう．
（ ここで T = 1000 付近では α = 7 程度となる. ）
しかしながら，より理論的には
x を適当な実数 ( π/6 ≧ x ≧ 0 ) として
sin(x), tan(x) のマクローリンの展開式で考えれば，
sin(x) = x - x3/3! + x5/5! … (4-3-1)
tan(x) = x + x3/3 + 2x5/15 … (4-3-2)
∴ tan(x) - sin(x) ≒ x3/2 (4-3-3)
この対数を取って
log(tan(x) - sin(x)) ≒ 3log(x) - log(2) (4-4-1)
を得る．
ここで x' = x/2 とすれば
log(tan(x') - sin(x')) ≒ 3log(x) - 4log(2) (4-4-2)
さらに 4-4-2 から 4-4-1 を引けば，
-3log(2) = -3*0.301…≒ -0.9 (4-4-3)
0.9 の前の - (マイナス) の意味は計算の上では小数点以下の精度の高い桁数の長い方への移動を意味していると考えられる．
すなわち，アルキメデスの方法で半角公式を一回実行する毎に
不等式
sin(x) ＜ x ＜ tan(x) (4-5)
の上限と下限が初めて不一致になる桁指数が約 0.9 桁づつより長くなって行くということだと考えられる．
これは半角公式一回実行直前の精度を基準にして考えると一回実行直後では約 0.9 桁精度が落ちることを意味していると言えるだろう．
言い換えると半角公式を実行する一段階前の桁数から既に約 0.9 桁ほど精度を高くしておかなければならなかったのだと考えられる．
この様に考えると一番最後の半角公式の実行から期間全体の初めまで遡れば，
3log(2)*T/(2log(2)) ≒ 3T/2 = 1.5T (4-6-1)
だけ桁数を長くして精度を上げておかなければならなかったことになる．
この結果，
最後の中心角における必要計算領域が約 1.5T 桁であったことも考慮に入れれば，
アルキメデスの方法における合理的な計算領域の総計は，
これを Sb とすれば，
Sb ≒ 1.5T + 1.5T = 3T (4-6-2)
であると考えられる．(一部誤り訂正と加筆 2013/4/20)
この Sb をアルキメデスの方法における理論的最小スケールと呼んでおこう．
（この桁数は実用的には数千桁付近だと約 30% 以上長すぎる。）
このような理由からアルキメデスの方法では後になるほど計算精度を低くして行くような意味の可変スケールは使えなくて，
初めから終わりまで固定精度で計算するしかないと考えられる．
アルキメデスの方法で要求精度が 1000 桁で固定スケールを 2007 桁とした場合
所要時間 0 分 41 秒 137 ミリ秒
アルキメデスの方法で要求精度が 1000 桁で固定スケールを 3000 桁とした場合
所要時間 1 分 43 秒 272 ミリ秒
この様にアルキメデスの方法では
半角公式を実行する毎に
前回の計算領域を今回の最小限度必要な計算領域よりも長くしておかなければ
計算の精度が落ちてしまうのである．
私はこれをアルキメデスのジレンマと呼んでいる．
これと比較して単位複素弧度回転演算子法でしかも二冪商の場合には
倍角公式を何回実行しても前回と同じ計算精度が得られるので計算の精度が落ちず，
しかも最小限度必要な計算領域も毎回ごとに短くなっていくので
計算領域を可変的に取ることができ計算速度も向上するのである．
（ これらに事がらで分かるように単位複素弧度回転演算子法でしかも二冪商の場合は
理論的な説明の合理性においても実用的な性能の面においても
アルキメデスの方法よりも優れていると考えられるのである．）
（５）アルキメデスの方法で平方根をニュートン法で計算した場合と単位複素弧度回転演算子法を応用した場合の
性能の比較とスケールの考察
アルキメデスの方法で要求精度が 1500 桁で平方根をニュートン法で計算した場合
所要時間 4 分 19 秒 537 ミリ秒
単位複素弧度回転演算子法で要求精度が 1500 桁の場合
所要時間 0 分 32 秒 885 ミリ秒
アルキメデスの方法で要求精度が 3000 桁で平方根をニュートン法で計算し経験的最小値の固定スケールとした場合
所要時間 17 分 28 秒 694 ミリ秒
単位複素弧度回転演算子法で要求精度が 3000 桁で要求精度の約 3/2 倍の桁数の固定スケールとした場合
所要時間 3 分 3 秒 877 ミリ秒
単位複素弧度回転演算子法で要求精度が 3000 桁で倍角公式 4 回毎に 1 桁づつ縮めるような可変スケールとした場合
所要時間 2 分 15 秒 782 ミリ秒
単位複素弧度回転演算子法で要求精度が 3000 桁で倍角公式 17 回毎に 5 桁づつ縮めるような可変スケールとした場合
所要時間 2 分 8 秒 325 ミリ秒
単位複素弧度回転演算子法で要求精度が 12000 桁で要求精度の約 3/2 倍の桁数の固定スケールとした場合
所要時間 3 時間 15 分 58 秒 298 ミリ秒
単位複素弧度回転演算子法で要求精度が 12000 桁で倍角公式 4 回毎に 1 桁づつ縮めるような可変スケールとした場合
所要時間 2 時間 19 分 47 秒 450 ミリ秒
単位複素弧度回転演算子法で要求精度が 12000 桁で倍角公式 17 回毎に 5 桁づつ縮めるような可変スケールとした場合
所要時間 2 時間 10 分 51 秒 761 ミリ秒
上記結果の検算として最後の繰り返しの分母 H に 1 を加算し求められた円周率を照合した場合
所要時間 2 時間 13 分 12 秒 363 ミリ秒
（ §２.２のみ 2011/08/05_Fr に追記 ）
§§３ 単位複素弧度回転演算子法と級数法の比較
§３.１ 公開研究課題
これをお読みの数学愛好家および研究者の皆様に以下のような公開課題を提出させて頂きます.
適当なプログラム言語を用いて以下の課題を実行し，結果を比較せよ.
（１）マチンの公式を用いて円周率を１０００桁近似し，実行時間を測定せよ．
クリックで実行例の文章ファイルが開けます
（２）オイラーの公式を用いて円周率を１０００桁近似し，実行時間を測定せよ．
クリックで実行例の文章ファイルが開けます
（３）単位複素弧度回転演算子法で二冪商と加法定理を応用して円周率を１０００桁近似し，実行時間を測定せよ．
クリックで実行例の文章ファイルが開けます
（４）単位複素弧度回転演算子法で三冪商と加法定理を応用して円周率を１０００桁近似し，実行時間を測定せよ．
クリックで実行例の文章ファイルが開けます
（５）さらに１００００桁で試して見よ．
マチンの公式による 10000 桁の近似
オイラーの公式による 10000 桁の近似
二冪商による 10000 桁の近似
三冪商による 10000 桁の近似
（６）グレゴリーの公式を用いて円周率をそれぞれ 7,8,9,10 桁近似し，実行時間を測定せよ．
グレゴリーの公式による 7 桁の近似
グレゴリーの公式による 8 桁の近似
グレゴリーの公式による 9 桁の近似
グレゴリーの公式による 10 桁の近似
（ ※ 些細なことだが PC で計算プログラムを実行する場合に電源プランによっても処理速度が変化するので注意を要する．
特にノート型 PC ではバッテリ駆動が原則なので，このことが特に当てはまると考えられる．
参考までに自作の円周率計算プログラム（ 円周率を 3083 桁計算し結果を文章ファイルにリダイレクトして出力する．）
を使って電源プラン毎の処理時間を測定した結果を示しておく．
[1] HP Optimized：バッテリー寿命（☆☆☆）パフォーマンス（☆☆☆）における所要時間
：222393[msec] = 3 分 42 秒 893 ミリ秒
[2] 省電力： バッテリー寿命（☆☆☆☆☆☆）パフォーマンス（☆☆）における所要時間
：577496[msec] = 9 分 37 秒 496 ミリ秒
[3] 高パーフォーマンス： バッテリー寿命（☆☆）パフォーマンス（☆☆☆☆☆☆） における所要時間
：221988[msec] = 3 分 41 秒 988 ミリ秒
３つの内では [3] が最も良いが，
しかし [1] と [3] の処理時間の比は 222393/221988 ≒ 1.001824 とほとんど同じである．
それで特に理由がない限り HP Optimized にしておけば良いと感じられる．
（ ただし HP Optimized は私個人の PC の呼称．他機種でも類似の機能があると思える．）
こんな訳で適当な長さの円周率の計算プログラムを自作しておいて自分の PC の性能を時たま調べて見るのも良いかも．^^）
（ 上記でオイラーの公式とマチンの公式を用いて 10000 桁の円周率を計算したときに
私が作成し使用した JaVa のソースプログラムのコピーを公開しておきます．(2011/5/7 追記)
オイラーの公式による 10000 桁の近似の Java ソースプログラム
マチンの公式による 10000 桁の近似の Java ソースプログラム
）
§３.２.１ 単位複素弧度回転演算子法の定義に必要な諸定理の初等的な立場からの証明
単位複素弧度回転演算子法の定義に必要な諸定理は我が国の高等学校卒業までに得られる事柄から証明可能である．
念のため補足事項としてこれらを証明しておこう．
（１）sin(θ) ＜ θ ＜ tan(θ) (1-1)
の証明
[図‐１]
与えられた任意の弧度を θ とせよ.
図‐１の如く原点 O を中心として半径 1, 中心角 θ で扇形 OAB を描け.
B から OA に垂線を下ろしその足を H とせよ.
A から OA 上に垂線 AT を立て，かつ，T は OB の延長線上にあるとせよ.
然らば各辺の長さは
OA = 1 (1-2-1)
HB = sin(θ) (1-2-2)
AT = tan(θ) (1-2-3)
ゆえに各部の面積において
三角形 OAB = (1/2)sin(θ) (1-3-1)
扇形 OAB = (1/2)θ (1-3-2)
三角形 OAT = (1/2)tan(θ) (1-3-3)
然るに
三角形 OAB ＜ 扇形 OAB ＜ 三角形 OAT (1-4-1)
∴ (1/2)sin(θ) ＜ (1/2)θ ＜ (1/2)tan(θ) (1-4-2)
∴ sin(θ) ＜ θ ＜ tan(θ)
//
（２）余弦と正弦の加法定理
cos(α±β) = cos(α)cos(β) -/+ sin(α)sin(β) (2-1-1)
sin(α±β) = sin(α)cos(β) ± cos(α)sin(β) (2-1-2)
の証明
[図‐２]
与えられたる任意の二つの弧度を α,β とせよ.
図‐２の如くに原点を O とし，O から任意の長さで水平に OA を引け．
A から OA に垂線 AB を立て，かつ，∠BOA = α とせよ.
B から OB に垂線 BC を立て, かつ, ∠COB = β とせよ.
C から OA に垂線 CE を下ろせ．
B から CE に垂線 BD を下ろせ．
然らば
∠ COA = α + β (2-2-1)
OB = OC・cos(β) (2-2-2)
BC = OC・sin(β) (2-2-3)
OA = OB・cos(α) = OC・cos(α)cos(β) (2-3-1)
AB = OB・sin(α) = OC・sin(α)cos(β) (2-3-2)
BD = BC・sin(α) = OC・sin(α)sin(β) (2-4-1)
CD = BC・cos(α) = OC・cos(α)sin(β) (2-4-2)
ゆえに
OE = OA - AE = OA - BD = OC・cos(α)cos(β) - OC・sin(α)sin(β) (2-5-1)
CE = CD + DE = CD + AB = OC・cos(α)sin(β) + OC・sin(α)cos(β) (2-5-2)
一方
OE = OC・cos(α+β) (2-6-1)
CE = OC・sin(α+β) (2-6-2)
ゆえに
cos(α+β) = cos(α)cos(β) - sin(α)sin(β)
sin(α+β) = sin(α)cos(β) + cos(α)sin(β)
さらに
cos(-β) = cos(β) (2-7-1)
sin(-β) = -sin(β) (2-7-2)
から
cos(α-β) = cos(α)cos(β) + sin(α)sin(β),
sin(α-β) = sin(α)cos(β) - cos(α)sin(β)
//
（３）ド・モアブルの定理
( cos(θ) + isin(θ))m = cos(mθ) + isin(mθ) (3-1)
の証明
(3-1) m = 1 で定理が真であるのは自明．
(3-2) m = n で定理が真であると仮定する．
すなわち，
( cos(θ) + isin(θ))n = cos(nθ) + isin(nθ) (3-2)
であることが既に示されているとせよ．
(3-3) m = n + 1 では，
定理の左辺は
( cos(θ) + isin(θ))n+1 = ( cos(θ) + isin(θ))n( cos(θ) + isin(θ)) (3-3-1)
ところが既に
( cos(θ) + isin(θ))n = cos(nθ) + isin(nθ)
が真であるので
( cos(θ) + isin(θ))n( cos(θ) + isin(θ)) = ( cos(nθ) + isin(nθ))( cos(θ) + isin(θ))
= ( cos(nθ)cos(θ) - sin(nθ)sin(θ)) + i( sin(nθ)cos(θ) + cos(nθ)sin(θ)) (3-3-2)
ところが加法定理によれば
cos(nθ)cos(θ) - sin(nθ)sin(θ) = cos((n+1)θ) (3-4-1)
sin(nθ)cos(θ) + cos(nθ)sin(θ) = sin((n+1)θ) (3-4-2)
∴ ( cos(θ) + isin(θ))n+1 = cos((n+1)θ) + isin((n+1)θ) (3-5)
ゆえに定理は m = n で真ならば m = n + 1 でも真である．
ところが定理は既に m = 1 で真であった．
ゆえに帰納法により定理は任意の m において真である．//
§３.２.２ 加法定理と複素数の乗法
x1,y1,x2,y2,x3,y3 を実数たち,
z1,z2,z3 を複素数たちとして,
z1 = x1 + iy1 (1-1)
z2 = x2 + iy2 (1-2)
z3 = x3 + iy3 (1-3)
z1・z2 = z3 (1-4)
であったとすれば,
z3 = (x1+iy1)(x2+iy2) = (x1x2-y1y2) + i(x1y2+y1x2)
∴ x3 = x1x2-y1y2 (2-1)
y3 = x1y2+y1x2 (2-2)
ところがここで適当な実数たち r1,r2,r3 弧度たちを θ1,θ2,θ3 として
いわゆる極座標形式でこれらを書き直してみれば,
r1 = √(x12+y12) (3-1)
r2 = √(x22+y22) (3-2)
r3 = √(x32+y32) (3-3)
θ1 = tan-1(y1/x1) (4-1)
θ2 = tan-1(y2/x2) (4-2)
θ3 = tan-1(y3/x3) (4-3)
x1 = r1cos(θ1) (5-1-1)
y1 = r1sin(θ1) (5-1-2)
x2 = r2cos(θ2) (5-2-1)
y2 = r2sin(θ2) (5-2-2)
x3 = r3cos(θ3) (5-3-1)
y3 = r3sin(θ3) (5-3-2)
z1 = r1(cos(θ1)+isin(θ1)) (6-1)
z2 = r2(cos(θ2)+isin(θ2)) (6-2)
z3 = r3(cos(θ3)+isin(θ3)) (6-3)
となることはよく知られている．
ゆえに, ここからさらに,
z3 = z1・z2
= r1(cos(θ1)+isin(θ1))・r2(cos(θ2)+isin(θ2))
= r1r2( ( cos(θ1)cos(θ2) - sin(θ1)sin(θ2)) + i( sin(θ1)cos(θ2) + cos(θ1)sin(θ2)) ) (6-4)
一方
cos(θ1)cos(θ2)-sin(θ1)sin(θ2) = cos(θ1+θ2) (7-1)
sin(θ1)cos(θ2)+cos(θ1)sin(θ2) = sin(θ1+θ2) (7-2)
などから,
r3 = r1r2 (8-1)
θ3 = θ1 + θ2 (8-2)
であることが分かる．
ここで特に,
r1 = r2 = 1
であったならば
Re(z1z2) = cos(θ1)cos(θ2)-sin(θ1)sin(θ2) = cos(θ1+θ2) (9-1)
Im(z1z2) = sin(θ1)cos(θ2)+cos(θ1)sin(θ2) = sin(θ1+θ2) (9-2)
となり，特に絶対値が 1 の複素数の乗法二項演算において，
実部の演算は余弦加法定理と同型であり, 虚部の演算は正弦加法定理と同型であると言える．
ゆえに，複素数がまだ定義されていない段階であっても, 余弦と正弦の加法定理の証明を与えて後，
虚数単位（ i ）と同様の規則だけを臨時に採用すれば単位複素弧度回転演算子法を定義することは可能であると言えよう．
§３.３.１ 適格な級数の有限部分和の下限と上限を用いた不等式により逆正接関数が近似できることの証明
級数法においては逆正接関数をテーラーの定理によって展開した公式が主に用いられている．
しかしながら収束の遅い級数の代表でもあり，また収束半径の真上で公式を適用しようとする級数の代表でもあるグレゴリーの公式
などにおいては公式が本当に成立するのだろうかという一抹の不安が常に付きまとう．
このような疑念が生じる理由は無限級数の収束を根拠として公式が等式で与えられていることもその一因であろうと思われる．
本節では無限や極限値の概念を表面的には避け有限の概念だけを用い，等式ではなく不等式を用いた公式の証明を掲げておこう．
（ しかしながら基本的な関数の微積分において無限や極限値の概念を用いることはやむを得ないことと考える．）
n を任意の自然数として次のような有限級数 s を考えよう．
s = ∑(1,n)(-1)n-1x2(n-1) = 1 - x2 + x4 - + … +(-1)n-1x2(n-1) (0＜x≦1) (1-1)
ここで最終項の形によって s を sa と sb の二種類に分けるものとする．
sa = ∑(1,n)(-1)n-1x2(n-1) (n=2m) (1-2)
sb = ∑(1,n)(-1)n-1x2(n-1) (n=2m+1) (1-3)
すると次の不等式が成立する．
◇ sa ＜ 1/(1+x2) ＜ sb (0＜x≦1) (1-4)
[証明]
(1+x2)sa = 1 + (-1)2m-1x4m ＜ 1 (1-5-1)
∴ sa ＜ 1/(1+x2) (0＜x≦1) (1-5-2)
(1+x2)sb = 1 + (-1)2mx4m+2 ＞ 1 (1-6-1)
∴ sb ＞ 1/(1+x2) (0＜x≦1) (1-6-2)
（証明終了）
ここで注意したいことは (1-4) は x=1 でも成立することである．（ グレゴリーの公式が成立することの根拠となる．）
さてここで s を x について項別に積分した関数 S を考えよう．すなわち，
S = ∫s dx = x - x3/3 + x5/5 -+ … +(-1)(n-1)x2n-1/(2n-1)
= ∑(1,n)(-1)n-1x2n-1/(2n-1) (0＜x≦1) (1-7)
同時にまた sa,sb も x について項別に積分した関数を考え順にSa,Sb としよう．すなわち，
Sa = ∑(1,n)(-1)n-1x2n-1/(2n-1) (n=2m) (1-8)
Sb = ∑(1,n)(-1)n-1x2n-1/(2n-1) (n=2m+1) (1-9)
ところが，ここで
∫1/(1+x2) dx = tan-1(x) (1-10)
を既知とすれば，
◇ Sa ＜ tan-1(x) ＜ Sb (0＜x≦1) (1-11)
が成立することが分かる．
(1-11) を以後の使用に便利なようにもう少し実用的に書き直しておけば，
x-x3/3+x5/5-+…(-1)2m-1x4m-1/(4m-1)＜tan-1(x)＜x-x3/3+x5/5-+…(-1)2mx4m+1/(4m+1) (1-11-1)
となる．
さらに以下の公式が成立する．
◇ Sa(m) ＜ Sa(m+1) ＜ tan-1(x) ＜ Sb(m+1) ＜ Sb(m) (1-12)
[証明]
Sa(m+1) - Sa(m) = (-1)2mx4m+1/(4m+1) + (-1)2m+1x4m+3/(4m+3)
= x4m+1/(4m+3) ((4m+3)/(4m+1) - x2) ＞ 0 (1-12-1)
Sb(m+1) - Sb(m) = (-1)2m+2x4m+5/(4m+5) + (-1)2m+1x4m+3/(4m+3)
= x4m+3/(4m+5) ( x2 - (4m+5)/(4m+3)) ＜ 0 (1-12-2)
（証明終了）
これをもう少し具体的に書き直せば，
Sa(1)＜Sa(2)＜…＜Sa(m)＜tan-1(x)＜Sb(m)＜…＜Sb(2)＜Sb(1) (1-12-3)
∴ lim(m→ω)S(m) = tan-1(x) (1-12-4)
このように考えることで従来のテーラーの定理を用いた結論とも一致する．
しかしながらグレゴリーの公式に関しては，特に実用的な立場から,
収束しない級数ではないか？と言う先述の結論とは逆の否定的な考え方も存在するように思われるのである．
以下でこれを示そう．
tan-1(x) の展開式を
lim(n→ω) x - x3/3 + x5/5 -+ … +(-1)n-1x2n-1/(2n-1) → tan-1(x) (1-13-1-1)
であるとし，この両辺を x で微分すれば，
lim(n→ω) 1 - x2 + x4 -+ … +(-1)n-1x2(n-1) → 1/(1+x2) (1-13-1-2)
となり，この (1-13-1-1),(1-13-1-2) で x = 1 を代入すれば，
1 - 1/3 + 1/5 -+ … +(-1)n-1/(2n-1) → π/4 (1-13-2-1)
1 - 1 + 1 -+ … +(-1)n-1 → 1/2 (?) (1-13-2-2)
となるが，この (1-13-2-2) は明らかに偽である．
(1-13-2-2) では下限の値は常に 0 で上限の値は常に 1 であって間に挟まれた 1/2 に限りなく近づいていくと考えるには
明らかに無理があり，下限と上限は殆どというよりは全く動かないのである．
（ 少なくともこのことから (1-13-2-1) の左辺と右辺が完全には同じ関数でないか，
さもなければ x = 1 では成立しないかのどちらかであるとだけは言えていると思われる．）
(1-13-2-1) では下限の値と上限の値が間に挟んだ π/4 の値に項の番号の増加と共に徐々に近づいて行くようにも考えられるが
実際に私の PC で計算させた限りにおいては超～遅いのである．
参考までにこの分野ではお馴染みの級数を用いた場合の要求精度と計算所要時間を掲げておけば，
オイラーの公式で要求精度 10000 桁に対する所要時間は約 13 分 41 秒 387 ミリ秒
マチンの公式で要求精度 10000 桁に対する所要時間は約 4 分 47 秒 133 ミリ秒
程度であったが，一方これらに対してグレゴリーの公式では
僅か要求精度 10 桁であっても約 4 時間 47 分 40 秒 27 ミリ秒 も掛かってしまい，
要求精度 20 桁においては私の概算が正しければ約 600 万年以上も掛かってしまうのである．
結局，グレゴリーの公式は実用的な意味では収束級数であるとは考えない方が良いのかも知れない．
（ テーラーやマクローリンの展開式で重要な前提が少なくとも二つある．それは，
（１）f(x) が x で何回でも微分できること．（微分連続性）
（２）lim(n→∞) xn → 0 を無条件で真とする．
tan-1(x) は確かに（１）を満たしはするが，グレゴリーの公式では x = 1 であるとするために（２）と矛盾する．
ただ，それでも x ≒ 1 - δ (δ→0+) と考えれば辛うじて収束級数であると見なせなくも無いとも考えられる．
いずれにしても微妙である． ）
§３.３.２ マチンの公式の証明
恐らく世界で一番簡単であると私が信じるマチンの公式の証明が以下である.
θ1, θ2 を弧度（ラジアン）として
θ1 = tan-1(1/5) (1-1)
θ2 = tan-1(1/239) (1-2)
と定義せよ.
すると
cos(θ1)+isin(θ1) = 5/√(26)+i/√(26) (2-1)
cos(θ2)+isin(θ2) = 239/√(2392+1)+i/√(2392+1) (2-2)
となる．
然らば，簡単な複素数の計算により
(5/√(26)+i/√(26))4 = (239/√(2392+1)+i/√(2392+1))(1/√2+i/√2)
= (119/169)+i(120/169) (3-1)
となることが確かめられる.
∴ 4θ1 - θ2 = π/4 (3-2)
ゆえに tan-1(1/5), tan-1(1/239) のテーラーの展開式からよく知られたマチンの公式を得る.//
（ ここでも複素弧度回転演算子を使った計算が冗長を省かせている.）
参考のためマチンの公式の級数展開も導いておこう．
tan-1(x) をテーラーの定理によって展開すれば，良く知られているように，
tan-1(x) = x - (1/3)x3 + (1/5)x5 - (1/7)x7 + - … (4-1)
となる．
これを，
tan-1(x) = ∑(-1)(n-1)/2(1/n)xn + Rn ( n = 1,3,5,… ) (4-2)
と書いておくことにする．（ Rn は剰余項 ）
ここでさらに A(m), B(m) をそれぞれ下限と上限の有理数として，
A(m) = ∑(-1)(n-1)/2(1/n)xn ( n = 1,3,…,4m-1 ) (4-3-1)
B(m) = A(m) + (1/(2m+1))x4m+1 (4-3-2)
であるとすれば，
A(m) ＜ tan-1(x) ＜ B(m) (4-3-3-1)
あるいは，
A(m) ＜ tan-1(x) ＜ A(m) + (1/(4m+1))x4m+1 (4-3-3-2)
となる．
それでは，これらの式たちから tan-1(1/5), tan-1(1/239) を実際に求めてみよう．
A1(m)=(1/5)-(1/3)(1/5)3+- … +(1/(n-1))(1/5)n-1-(1/n)(1/5)n (n=4m-1) (4-4-1)
B1(m) = A1(m) + (1/(4m+1))(1/5)4m+1 (4-4-2)
となるので，
A1(m) ＜ tan-1(1/5) ＜ B1(m) (4-4-3-1)
あるいは，
A1(m) ＜ tan-1(1/5) ＜ A1(m) + (1/(4m+1))(1/5)4m+1 (4-4-3-2)
同様に
A2(m')=(1/239)-(1/3)(1/239)3+- … +(1/(n'-1))(1/239)n'-1-(1/n')(1/239)n' (n'=4m'-1) (4-5-1)
B2(m') = A2(m') + (1/(4m'+1))(1/239)4m'+1 (4-5-2)
から
A2(m') ＜ tan-1(1/239) ＜ B2(m') (4-5-3-1)
あるいは，
A2(m') ＜ tan-1(1/239) ＜ A2(m') + (1/(4m'+1))(1/239)4m'+1 (4-5-3-2)
となる．
ゆえに
4A1(m) ＜ 4tan-1(1/5) ＜ 4A1(m) + (4/(4m+1))(1/5)4m+1 (4-6-1)
-A2(m') - (1/(4m'+1))(1/239)4m'+1 ＜ -tan-1(1/239) ＜ -A2(m') (4-6-2)
4tan-1(1/5) - tan-1(1/239) = π/4 (4-6-3)
から，
4A1(m)-A2(m')-(1/(4m'+1))(1/239)4m'+1 ＜ π/4 ＜ 4A1(m)+(4/(4m+1))(1/5)4m+1-A2(m') (4-6-4)
が得られる．そこで新たに，
C(m,m') = 8A1(m)-2A2(m')-(2/(4m'+1))(1/239)4m'+1 (4-7-1)
D(m,m') = 8A1(m)-2A2(m')+(8/(4m+1))(1/5)4m+1 (4-7-2)
とすれば，
C(m,m') ＜ π/2 ＜ D(m,m') (4-7-3)
となるので，π の近似値を Pi(m,m'), 誤差を Er(m,m') とすれば，
Pi(m,m') = D(m,m') + C(m,m') (4-8-1)
Er(m,m') = D(m,m') - C(m,m') = (8/(4m+1))(1/5)4m+1 + (2/(4m'+1))(1/239)4m'+1 (4-8-2)
からマチンの公式による円周率の近似が得られる．
ここで注意すべきことは，近似の要求精度が T であったときにもはや m と m' の最小値は等しくないということである．
この m,m' を求めておこう．
このとき精度が T であるような tan-1(1/5) の級数展開において，
4B1(m) - 4A1(m) = (4/(4m+1))(1/5)4m+1 (4-9-1)
log10((4/(4m+1))(1/5)4m+1) ≒ -T (4-9-2)
から
(4m)log105 ≒ T (4-9-3)
∴ m ≒ T/(4log105) ≒ T/(4・0.699…) ≒ 0.3577T (4-9-4)
およその項の個数を N1 とすれば，
N1 ≒ 2m ≒ 0.7153T (4-9-5)
同様に精度が T であるような tan-1(1/239) の級数展開において，
B2(m') - A2(m') = (1/(4m'+1))(1/239)4m'+1 (4-9-1)
log10((1/(4m'+1))(1/239)4m'+1) ≒ -T (4-9-2)
から
(4m')log10239 ≒ T (4-9-3)
∴ m' ≒ T/(4log10239) ≒ T/(4・2.378…) ≒ 0.1051T (4-9-4)
およその項の個数を N2 とすれば，
N2 ≒ 2m' ≒ 0.2102T (4-9-5)
それで近似式全体においてのおよその項の個数を Ns とすれば，
Ns ≒ N1 + N2 ≒ 0.7153T + 0.2102T = 0.9255T (4-9-5)
となる．
一方，これと比較して単位複素弧度回転演算子法で二冪商のときの複素二項演算のおよその回数 Nc は
Nc ≒ T/(2log102) ≒ 1.661T (4-10-1)
となる．
この結果，前者の項の個数 Ns と後者の複素二項演算の回数 Nc との単純な比は，
Ns/Nc ≒ 0.9255/1.661 ≒ 0.55719 (4-10-2)
となる．
私の最近の実行例では要求精度 10000 桁に対して
前者の場合の所要時間 hs = 287133 [msec] = 4 分 47 秒 133 ミリ秒
後者の場合の所要時間 hc = 5255656 [msec] = 1 時間 27 分 35 秒 656 ミリ秒
所要時間の比 hs/hc ≒ 287133/5255656 ≒ 0.0546331
これは同じ要求精度に対する所要時間の意味で前者は後者の 約 1/20 となることを示していると考えられる．
前者の項一回分の平均処理時間 hs/Ns ≒ 287133/(0.9255・10000) ≒ 31.0246 [msec/回]
後者の複素数乗法一回分の平均処理時間 hc/Nc ≒ 5255656/(1.661・10000) ≒ 316.415 [msec/回]
一回あたりの処理時間の比 31.0246/316.415 ≒ 0.098050
この数値は平均的な処理時間の意味で前者の一回分は後者の一回分の約 1/10 であることを示していると考えられる．
§３.３.３ （ 級数法による円周率の計算における ）オイラーの公式の証明
オイラーの公式もマチンの公式と原理的には類似している. これを示そう．
θ1, θ2 を弧度（ラジアン）として
θ1 = tan-1(1/2) (1-1)
θ2 = tan-1(1/3) (1-2)
と定義せよ.
すると
cos(θ1)+isin(θ1) = 2/√5+i/√5 (2-1)
cos(θ2)+isin(θ2) = 3/√(10)+i/√(10) (2-2)
となる.
然らば，簡単な複素数の計算により
(2/√5+i/√5)(3/√(10)+i/√(10)) = (1/√2)+i(1/√2) (3-1)
となることが確かめられる.
∴ θ1 + θ2 = π/4 (3-2)
すなわち
tan-1(1/2) + tan-1(1/3) = π/4 (3-3)
ゆえに tan-1(1/2), tan-1(1/3) のテーラーの展開式からオイラーの公式を得る.//
参考のためオイラーの公式の級数展開も導いておこう．
tan-1(x) をテーラーの定理によって展開すれば，良く知られているように，
tan-1(x) = x - (1/3)x3 + (1/5)x5 - (1/7)x7 + - … (4-1)
となる．
これを
tan-1(x) = ∑(-1)(n-1)/2(1/n)xn + Rn ( n = 1,3,5,… ) (4-2)
と書いておくことにする．（ Rn は剰余項 ）
ここでさらに A(m), B(m) をそれぞれ下限と上限の有理数として，
A(m) = ∑(-1)(n-1)/2(1/n)xn ( n = 1,3,…,4m-1 ) (4-3-1)
B(m) = A(m) + (1/(2m+1))x4m+1 (4-3-2)
であるとすれば，
A(m) ＜ tan-1(x) ＜ B(m) (4-3-3-1)
あるいは，
A(m) ＜ tan-1(x) ＜ A(m) + (1/(4m+1))x4m+1 (4-3-3-2)
となる．
それでは，これらの式たちから tan-1(1/2), tan-1(1/3) を実際に求めてみよう．
A1(m)=(1/2)-(1/3)(1/2)3+- … +(1/(n-1))(1/2)n-1-(1/n)(1/2)n (n=4m-1) (4-4-1)
B1(m) = A1(m) + (1/(4m+1))(1/2)4m+1 (4-4-2)
となるので，
A1(m) ＜ tan-1(1/2) ＜ B1(m) (4-4-3-1)
あるいは，
A1(m) ＜ tan-1(1/2) ＜ A1(m) + (1/(4m+1))(1/2)4m+1 (4-4-3-2)
同様に
A2(m')=(1/3)-(1/3)(1/3)3+- … +(1/(n'-1))(1/3)n'-1-(1/n')(1/3)n' (n'=4m'-1) (4-5-1)
B2(m') = A2(m') + (1/(4m'+1))(1/3)4m'+1 (4-5-2)
となるので，
A2(m') ＜ tan-1(1/3) ＜ B2(m') (4-5-3-1)
あるいは，
A2(m') ＜ tan-1(1/3) ＜ A2(m') + (1/(4m'+1))(1/3)4m'+1 (4-5-3-2)
となる．さらに，
A1(m) ＜ tan-1(1/2) ＜ A1(m) + (1/(4m+1))(1/2)4m+1 (4-6-1)
A2(m') ＜ tan-1(1/3) ＜ A2(m') + (1/(4m'+1))(1/3)4m'+1 (4-6-2)
tan-1(1/2) + tan-1(1/3) = π/4 (4-6-3)
から，
A1(m)+A2(m') ＜ π/4 ＜ A1(m)+A2(m')+(1/(4m+1))(1/2)4m+1+(1/(4m'+1))(1/3)4m'+1 (4-6-4)
が得られる．そこで新たに，
C(m,m') = 2A1(m)+2A2(m') (4-7-1)
D(m,m') = 2A1(m)+2A2(m')+(2/(4m+1))(1/2)4m+1+(2/(4m'+1))(1/3)4m'+1 (4-7-2)
とすれば，
C(m,m') ＜ π/2 ＜ D(m,m') (4-7-3)
となるので，π の近似値を Pi(m,m'), 誤差を Er(m,m') とすれば，
Pi(m,m') = D(m,m') + C(m,m') (4-8-1)
Er(m,m') = D(m,m') - C(m,m') = (2/(4m+1))(1/2)4m+1 + (2/(4m'+1))(1/3)4m'+1 (4-8-2)
からオイラーの公式による円周率の近似が得られる．
ここでも注意すべきことは，近似の要求精度が T であったときにもはや m と m' の最小値は等しくないということである．
この m,m' を求めておこう．
このとき精度が T であるような tan-1(1/2) の級数展開において
B1(m) - A1(m) = (1/(4m+1))(1/2)4m+1 (4-9-1)
log10((1/(4m+1))(1/2)4m+1) ≒ -T (4-9-2)
から
(4m)log102 ≒ T (4-9-3)
∴ m ≒ T/(4log102) ≒ T/(4・0.301…) ≒ 0.8306T (4-9-4)
およその項の個数を N1 とすれば
N1 ≒ 2m ≒ 1.6612T (4-9-5)
同様に精度が T であるような tan-1(1/3) の級数展開において
B2(m') - A2(m') = (1/(4m'+1))(1/3)4m'+1 (4-9-1)
log10((1/(4m'+1)(1/3)4m'+1) ≒ -T (4-9-2)
から
(4m')log103 ≒ T (4-9-3)
∴ m' ≒ T/(4log103) ≒ T/(4・0.477…) ≒ 0.5240T (4-9-4)
およその項の個数を N2 とすれば
N2 ≒ 2m' ≒ 1.0480T (4-9-5)
それで近似式全体においてのおよその項の個数を Ns とすれば
Ns ≒ N1 + N2 ≒ 1.6612T + 1.0480T = 2.709T (4-9-5)
となる．
一方，これと比較して単位複素弧度回転演算子法で二冪商のときの複素二項演算のおよその回数 Nc は
Nc ≒ T/(2log102) ≒ 1.661T (4-10-1)
となる．
この結果，前者の項の個数 Ns と後者の複素二項演算の回数 Nc との単純な比は，
Ns/Nc ≒ 2.709/1.661 ≒ 1.631 (4-10-2)
となる．
私の最近の実行例では要求精度 10000 桁に対して
前者の場合の所要時間 hs = 821387 [msec] = 13 分 41 秒 387 ミリ秒
後者の場合の所要時間 hc = 5255656 [msec] = 1 時間 27 分 35 秒 656 ミリ秒
所要時間の比 hs/hc ≒ 821387/5255656 ≒ 0.156286
これは同じ要求精度に対する所要時間の意味で前者は後者の 約 1/6 となることを示していると考えられる．
前者の項一回分の平均処理時間 hs/Ns ≒ 821387/(2.709・10000) ≒ 30.3206 [msec/回]
後者の複素数乗法一回分の平均処理時間 hc/Nc ≒ 5255656/(1.661・10000) ≒ 316.4 [msec/回]
一回あたりの処理時間の比 30.3206/316.4 ≒ 0.09582996
この数値は平均的な処理時間の意味で前者の一回分は後者の一回分の約 1/10 であることを示していると考えられる．
また 10000 桁の要求精度の下でマチンの方法とオイラーの方法との比較は
項数の マチン/オイラー ≒ 9255/27090 ≒ 0.34163898 ≒ 1/3
所要時間の マチン/オイラー = 287133/821387 ≒ 0.34957 ≒ 1/3
となっている．
§３.３.４ オイラーの方法やマチンの方法の一般化
オイラーの方法は実は容易に一般化でき，従って無数に多くの類似な方法で円周率の近似が可能である．
次の公式が成り立つ．
[オイラーの方法の一般原理]
a,b を任意の互いに素な二つの整数とせよ．
一般に次の公式が成り立つ．
tan-1(b/a)+tan-1((a-b)/(a+b))=π/4 (1)
[証明]
tan-1(b/a)=θ1 (2-1)
tan-1((a-b)/(a+b))=θ2 (2-2)
とせよ．然らば，
cos(θ1)+isin(θ1)=(a/√(a2+b2))+i(b/√(a2+b2)) (3-1)
cos(θ2)+isin(θ2)=((a+b)/√((a+b)2+(a-b)2))+i((a-b)/√((a+b)2+(a-b)2))
=((a+b)/√(2(a2+b2))+i((a-b)/√(2(a2+b2)) (3-2)
∴ (cos(θ1)+isin(θ1))(cos(θ2)+isin(θ2))
= ((a/√(a2+b2))+i(b/√(a2+b2)))(((a+b)/√(2(a2+b2))+i((a-b)/√(2(a2+b2)))
= ((a(a+b)-b(a-b))+i(a(a-b)+b(a+b)))/((√2)(a2+b2))
= ((a2+b2)+i(a2+b2))/((√2)(a2+b2))
= (1/√2)+i(1/√2) (4)
∴ θ1+θ2=π/4 (5)
∴ tan-1(b/a)+tan-1((a-b)/(a+b))=π/4 //
例えば, ここで a=2, b=1 とすれば
a+b=2+1=3
a-b=2-1=1
となるので
tan-1(1/2)+tan-1(1/3)=π/4
が成立し，オイラーの方法となる．
ここで要求精度を T (桁)とした場合, テーラー展開による近似計算に必要な項の個数 N はおよそ
N ≒ (T/2)(1/log10(a/b)+1/log10(a'/b'))
= (T/2)(1/(log10a-log10b)+1/(log10a'-log10b')) (6-1)
となる．
∴ N ∝ 1/(log10a-log10b)+1/(log10a'-log10b') (6-2)
それで
1/(log10a-log10b)+1/(log10a'-log10b') (6-3)
の値が小さいほど項の個数 N が少なくなると考えられる．
結局，
|b/a|, |b'/a'| が小さいほど項の個数が少なくなると考えられる．
理想的には |b|=|b'|=1 が望ましい．
b=1 は自然に満たせるが，普通の方法では｜b'|=1 にすることは不可能である．
ところが, オイラーの方法ではそれが満たされているのである．
オイラーの方法の場合
1/(log10(2/1))+1/(log10(3/1))) = 1/(log102)+1/(log103)
≒ 1/0.3+1/0.5 ≒ 5.3 (7-1)
となっている．
これがマチンの方法では
1/(log10(5/1))+1/(log10(239/1)) = 1/(log105)+1/(log10239)
≒ 1/0.7+1/2.37 ≒ 1.8 (7-2)
とオイラーの約 1/3 倍ほどとなってさらに理想的な値となっているのである．
（ 私の最近の実行結果では 10000 桁の要求精度に対して
オイラーの方法では約 13 分 42 秒未満，マチンの方法では約 4 分 48 秒未満であった．）
マチンの方法においては
a=5,b=1 (8-1)
でしかも
4tan-1(b/a)=tan-1(v/u) (8-2)
とすれば，
v/u = 120/119 (8-3)
a'=u+v=119+120=239 (8-4)
b'=u-v=119-120=-1 (8-5)
のようになっているのである．
オイラーにしてもマチンにしても稀な理想を見事に見つけているのであって，その眼力は天才的である．
マチンの公式による 10000 桁の近似
オイラーの公式による 10000 桁の近似
適当に a,b,u,v,a',b' を定めて代用の級数を見つける場合に，マチンより早く収束する級数を見つけることは難しいけれども，
オイラーと比べてほんの少しより良い級数を見つけることは可能である
例えば，私が自ら手計算によって発見した代用品として
2tan-1(1/3) + tan-1(1/7) = π/4 (9-1)
があったが，残念ながらこれは ハットンの公式と呼ばれて既に発見されていた．
ハットンの公式 (9-1) による 10000 桁の近似
ハットンはまた
3tan-1(1/4) + tan-1(5/99) = π/4 (9-2)
のような公式も発見している．
ハットンの公式 (9-2)による 10000 桁の近似
分子が１を超えても良い，分数に使用される整数の桁数が巨大であっても良いというのなら
私手作りの検算プログラムで発見した以下のような代用級数がある．
いづれも要求精度 10000 桁の円周率の計算を 8 分以内に完了した．
（ これらはいづれもオイラーを凌いではいるもののマチンには及ばなかった．）
(1) 発見一号
8tan-1(7/71) + tan-1((U-V)/(U+V)) = π/4 (10-1-1)
U+V=949261769468944, (10-1-2)
U-V=-753203917808, (10-1-3)
(U+V)/(U-V)≒1260.29850, (10-1-4)
発見一号による 10000 桁の近似
(2) 発見二号
16tan-1(3/61) + tan-1((U-V)/(U+V)) = π/4 (10-2-1)
U+V=52988971989839953101362774272, (10-2-2)
U-V=-45230722336591359678922496, (10-2-3)
(U+V)/(U-V)≒-1171.52610 (10-2-4)
発見二号による 10000 桁の近似
(3) 発見三号
8tan-1(111/1127) + tan-1((U-V)/(U+V)) = π/4 (10-3-1)
U+V=3825399223814540323147792, (10-3-2)
U-V=-5067212416336591856, (10-3-3)
(U+V)/(U-V)≒-754931.688176 (10-3-4)
発見三号による 10000 桁の近似
（ PC を利用して特に，級数法で円周率を求める場合，
円周率を計算することよりも，むしろ優れた級数を発見することにこそ PC を利用すべきなのかも知れない．）
最近の私の発見としてマチンの公式
4tan-1(1/5) - tan-1(1/239) = π/4 (11-1)
において，
tan-1(1/5) = tan-1(1/7) + tan-1(1/18) (11-2)
で式を書き改めると
4tan-1(1/7) + 4tan-1(1/18) - tan-1(1/239) = π/4 (11-3)
となるという発見があった．
マチン改一号による 10000 桁の近似
これだと 1/7 や 1/18 が 1/5 より小さい分数なので従来のマチンの公式よりは早めに収束するのではと期待していたが，
結果はむしろ若干遅めになって改善には至らなかった．
また
tan-1(1/5) = tan-1(1/6) + tan-1(1/31) (11-4)
のような置き換えも可能である．
この式を用いてもマチンの公式を変形することはできる．
マチン改ニ号による 10000 桁の近似
これも先ほどのマチン改一号とほぼ同程度の性能であって従来のマチンの公式には及ばなかった．
ここでの置き換えの一般原理は次のようになる．
一般に
tan-1(1/5) = tan-1(c/b) + tan-1(1/a) (11-5-1)
とするには n を適当な自然数として
a = 13n + 5 (11-5-2)
b = Re((5+i)(a-i))/k (11-5-3)
c = Im((5+i)(a-i))/k (11-5-4)
k = (kb,kc) (11-5-5)
とすれば良い．
同様に tan-1(1/239) についても
tan-1(1/239) = tan-1(c/b) + tan-1(1/a) (11-6-1)
とするには n を適当な自然数として
a = 169n + 239 (11-6-2)
b = Re((239+i)(a-i))/k (11-6-3)
c = Im((239+i)(a-i))/k (11-6-4)
k = (kb,kc) (11-6-5)
とすれば良い．
例えば，n=1 なら
tan-1(1/239) = tan-1(1/577) + tan-1(1/408) (11-6-6)
が成り立ち，これを先述の 11-3 に代入すると
4tan-1(1/7) + 4tan-1(1/18) - tan-1(1/408) - tan-1(1/577) = π/4 (11-6-6)
が得られる．
マチン改三号による 10000 桁の近似
これも性能においては従来のマチンの公式に及ばなかった．
さらに別の置き換えとして
tan-1(1/7) = tan-1(1/9) + tan-1(1/32) (11-6-6)
が成り立つので，これを先述の 11-3 に代入すると
4tan-1(1/9) + 4tan-1(1/18) + 4tan-1(1/32) - tan-1(1/239) = π/4 (11-6-7)
が得られる．
マチン改四号による 10000 桁の近似
これも性能においては従来のマチンの公式に及ばなかった．
が，しかし複数の PC で異なる tan-1 部を計算して，後で集計するというような方法でも採用すれば
改良されたマチンの公式のどれもが従来のものより若干速くなるだろうと考えられる．
さらに，
tan-1(1/9) = tan-1(1/11) + tan-1(1/50) (11-6-8)
から
4tan-1(1/11) + 4tan-1(1/18) + 4tan-1(1/32) + 4tan-1(1/50)- tan-1(1/239) = π/4 (11-6-9)
が得られる．
マチン改五号による 10000 桁の近似
これも性能においては従来のマチンの公式に及ばなかったが，複数の PC で同時に計算して後で集計するような場合には
従来のマチンの公式より速くなるはずである．
なぜなら tan-1 の部分で使用される最大の分数が従来のマチンの公式では 1/5 だが (11-6-9) では 1/11 だからである．
ここまでの例から経験的に言えることは
マチンの公式は単独の PC で円周率を近似計算する場合の最も簡単で最も速い最適解であろうということである．
§３.３.５ 逆正接関数の加法推移分解によるマチンの公式の変形と連弾法による円周率計算の高速化
次の公式が成り立つ．
◇ a を自然数かつ a ≧ 2 として
tan-1(2/a2) = tan-1(1/(a-1)) - tan-1(1/(a+1)) (1-1)
= tan-1(1/(a2-a+1)) + tan-1(1/(a2+a+1)) (1-2)
= tan-1(1/((a-1)2+(a-1)+1)) + tan-1(1/((a+1)2-(a+1)+1)) (1-3)
[証明]
tan-1(1/(a-1)) - tan-1(1/(a+1))
⇒ ((a-1)+i)((a+1)-i) = (a2-1+1)+i((a+1)-(a-1)) = a2+2i
∴ (1-1)が証明された．
tan-1(1/(a2-a+1)) + tan-1(1/(a2+a+1))
⇒ ((a2-a+1)+i)((a2+a+1)+i)
= a4+a2+2(a2+1)i = (a2+1)(a2+2i)
∴ (1-2)が証明された．
a2-a+1=(a-1)2+(a-1)+1
∴ (1-3)が証明された．//
例えばこれをマチンの公式の tan-1(1/5) に応用したならば，
tan-1(1/5)
=tan-1(1/7)+tan-1(1/18)
=tan-1(1/9)+tan-1(1/18)+tan-1(1/32)
=tan-1(1/11)+tan-1(1/18)+tan-1(1/32)+tan-1(1/50)
=tan-1(1/13)+tan-1(1/18)+tan-1(1/32)+tan-1(1/50)+tan-1(1/72)
・ ・ ・
=tan-1(1/25)+tan-1(1/18)+ … +tan-1(1/(24*12))
・ ・ ・
=tan-1(1/(5+2n))+tan-1(2/(4+2)2)+ … +tan-1(2/(4+2n)2) (n=1,2,3,…) (2)
が際限なく成立する．
あるいは更に巧妙な応用としてハットンの公式を変形して後，マチンの公式を当てはめると次のような公式が成立する．
◇ 6tan-1(1/8) + 2tan-1(1/57) + tan-1(1/239) = π/4 (3-1)
[証明]
ハットンの第一公式から
2tan-1(1/3) + tan-1(1/7) = tan-1(1) (3-2)
ところが
tan-1(1/3) = tan-1(1/5) + tan-1(1/8) (3-3)
( ∵ (3+i)(5-i) = 2(8+i) )
∴ 2tan-1(1/5) + 2tan-1(1/8) + tan-1(1/7) = tan-1(1) (3-4)
ところが
tan-1(1/7) = tan-1(1/8) + tan-1(1/57) (3-3)
( ∵ (7+i)(8-i) = 57+i, 57 = 72+7+1 )
∴ 2tan-1(1/5) + 3tan-1(1/8) + tan-1(1/57) = tan-1(1)
さらに両辺を二倍して
4tan-1(1/5) + 6tan-1(1/8) + 2tan-1(1/57) = 2tan-1(1) (3-5)
ところがマチンの公式
4tan-1(1/5) - tan-1(1/239) = tan-1(1) (3-6)
から
6tan-1(1/8) + 2tan-1(1/57) + tan-1(1/239) = π/4
（証明終了）
この公式 (3-1) は既に公に知られているようである．（ 証明だけは私独自のものであるが．）
さて, ここで PC が二台あってしかも社内 LAN のようにネットワーク化されていると仮定しよう．
そこで例えばマチンの公式により級数法で円周率を計算する場合に次のような方法を用いたとする．
(0) 計算の要求精度を T（桁）とする．
(1) 一台目の PC で tan-1(1/239) を除いた部分を実行する．
(2) 二台目の PC で tan-1(1/5) を除いた部分を実行する．
(3) 両者の処理のいづれか遅い方が終了した時点で終結不等式を導き円周率の計算を完了する．
こうすることで共通な前後の補助的処理と主として tan-1(1/5) の処理に要した時間程で計算を完了するはずである．
これを一台の PC だけを用いて測定するのに私は以下のような方法を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/239) を除いた部分を実行する．
(3) flg==2 で tan-1(1/5) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 285464[ミリ秒] = 4 分 45 秒 464 ミリ秒
(2) flg==1 のとき 220834[ミリ秒] = 3 分 40 秒 834 ミリ秒
(3) flg==2 のとき 65566[ミリ秒] = 1 分 5 秒 566 ミリ秒
同様に上の公式 (3-1) においても次のような方法を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/57) と tan-1(1/239) を除いた部分を実行する．
(3) flg==2 で tan-1(1/8) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 318583[ミリ秒] = 5 分 18 秒 583 ミリ秒
(2) flg==1 のとき 166640[ミリ秒] = 2 分 46 秒 640 ミリ秒
(3) flg==2 のとき 153270[ミリ秒] = 2 分 33 秒 270 ミリ秒
どうやら二台の PC で手分けして計算させた場合はマチンの公式より公式 (3-1) の方が約 1 分ほど速くなるようである．
また公式 (3-1) では二台の PC が同程度の負担となるようである．
さらに次の公式を私流の計算法から導いた．
◇ 12tan-1(1/18) + 8tan-1(1/57) - 5tan-1(1/239) = π/4 (4-1)
（ この公式はガウスによって既に発見されていた．）
[証明]
tan-1(1/5) = tan-1(1/7) + tan-1(1/18) (4-2)
( ∵ (5+i)(7-i)=2(18+i) )
tan-1(1/7) = tan-1(1/8) + tan-1(1/57) (4-3)
( ∵ (7+i)(8-i)=57+i) )
∴ tan-1(1/5) = tan-1(1/8) + tan-1(1/18) + tan-1(1/57)
∴ 4tan-1(1/5) = 4tan-1(1/8) + 4tan-1(1/18) + 4tan-1(1/57) (4-4)
ところがマチンの公式から
4tan-1(1/5) - tan-1(1/239) = tan-1(1) (4-5)
∴ 4tan-1(1/8) + 4tan-1(1/18) + 4tan-1(1/57) - tan-1(1/239) = tan-1(1) (4-6)
先述の (3-1) から
6tan-1(1/8) + 2tan-1(1/57) + tan-1(1/239) = tan-1(1) (3-1)
ゆえに (4-6) の両辺の３倍から (3-1) の両辺の２倍を引くことで
12tan-1(1/18) + 8tan-1(1/57) - 5tan-1(1/239) = π/4
（証明終了）
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/57) と tan-1(1/239) を除いた部分を実行する．
(3) flg==2 で tan-1(1/18) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 274716[ミリ秒] = 4 分 34 秒 716 ミリ秒
(2) flg==1 のとき 122709[ミリ秒] = 2 分 2 秒 709 ミリ秒
(3) flg==2 のとき 153036[ミリ秒] = 2 分 33 秒 36 ミリ秒
これは驚くべきことに一台の PC で実行した場合であってもマチンの公式より約 10 秒ほど速い．
従って今まで述べた級数の中ではどれよりも速い．
マチンの公式による 10000 桁の近似 (flg==0)
(3-1) による 10000 桁の近似 (flg==0)
(4-1) による 10000 桁の近似 (flg==0)
さらに次の公式も成り立つ．
◇ 20tan-1(1/57) + 12tan-1(3/79) - 5tan-1(1/239) = π/4 (5-1)
[証明]
先述の (4-1) から
12tan-1(1/18) + 8tan-1(1/57) - 5tan-1(1/239) = π/4 (4-1)
ところが
tan-1(1/18) = tan-1(1/57) + tan-1(3/79) (4-2)
( ∵ (18+i)(57-i) = 13(79+3i) )
∴ 20tan-1(1/57) + 12tan-1(3/79) - 5tan-1(1/239) = π/4
（証明終了）
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(3/79) と tan-1(1/239) を除いた部分を実行する．
(3) flg==2 で tan-1(1/57) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 277399[ミリ秒] = 4 分 37 秒 399 ミリ秒
(2) flg==1 のとき 88390[ミリ秒] = 1 分 28 秒 390 ミリ秒
(3) flg==2 のとき 190164[ミリ秒] = 3 分 10 秒 164 ミリ秒
これは一台の PC で実行した場合に残念ながら (4-1) より 2 秒ほど遅いがマチンの公式よりは約 8 秒ほど速い．
(4-1) より若干遅れる理由として tan-1(3/79) の分子の 3 の冪の処理が案外長く掛かるためと考えられる．
必ずしも性能の改善になるとは限らないが (5-1) を変形することはできる．
次の定理が成り立つ．（ この公式を私は正規化公式と呼んでいる．）
◇ tan-1(b/(ab±c)) = tan-1(1/a)-(±)tan-1(c/(a(ab±c)+b)) (6-1)
( ∵ (ab±c+ib)(a-i) = (a(ab±c)+b-i(ab-ab-(±)c = a(ab±c)+b -(±)ic )
ここで (a,b,±c) = (26,3,-1) として (6-1) を (5-1) に代入すると
◇ 12tan-1(1/26) + 20tan-1(1/57) - 5tan-1(1/239) - 12tan-1(1/2057) = π/4 (7-1)
を得る．( 3*262 + 26 + 3 = 2057 )
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/57) と tan-1(1/239) と tan-1(1/2057)を除いた部分を実行する．
(3) flg==2 で tan-1(1/26) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 307320[ミリ秒] = 5 分 7 秒 320 ミリ秒
(2) flg==1 のとき 108701[ミリ秒] = 1 分 48 秒 701 ミリ秒
(3) flg==2 のとき 199805[ミリ秒] = 3 分 19 秒 805 ミリ秒
これは一台の PC で実行した場合に残念ながらマチンの公式よりも約 25 秒ほども遅いので改善にはならなかった．
しかし二台以上の PC で実行する場合には良いかも知れない．式の形も綺麗である．
さらに
(57+i)(239-i)=26(524+i7)
から
◇ 12tan-1(3/79) + 20tan-1(7/524) + 15tan-1(1/239) = π/4 (8-1)
を得る．
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(7/524) と tan-1(1/239) を除いた部分を実行する．
(3) flg==2 で tan-1(3/79) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 287617[ミリ秒] = 4 分 47 秒 617 ミリ秒
(2) flg==1 のとき 125408[ミリ秒] = 2 分 5 秒 408 ミリ秒
(3) flg==2 のとき 163706[ミリ秒] = 2 分 43 秒 706 ミリ秒
これは一台の PC で実行した場合にマチンの公式よりも約 2 秒ほど遅いだけで殆ど対等であると言えよう．
さらに次の公式も私が発見した．
◇ 16tan-1(1/21) + 3tan-1(1/239) + 4tan-1(3/1042) = π/4 (9-1)
[証明]
ハットンの第二公式から
3tan-1(1/4) + tan-1(5/99) = tan-1(1) (9-2)
ところが
tan-1(1/4) = tan-1(1/5) + tan-1(1/21) (9-3)
( ∵ (4+i)(5-i) = 21+i )
tan-1(5/99) = tan-1(1/20) + tan-1(1/1985) (9-4)
（ 先述の正規化公式 (6-1) で (a,b,±c）= (20,5,-1) を代入 ）
∴ 12tan-1(1/5) + 12tan-1(1/21) + 4tan-1(1/20) + 4tan-1(1/1985) = 4tan-1(1) (9-5)
ところがマチンの公式から
12tan-1(1/5) - 3tan-1(1/239) = 3tan-1(1) (9-6)
(9-5) から (9-6) を辺々差し引けば
12tan-1(1/21) + 4tan-1(1/20) + 4tan-1(1/1985) + 3tan-1(1/239) = tan-1(1) (9-7)
ところが
tan-1(1/20) - tan-1(1/21) = tan-1(1/421) (9-8)
( ∵ (20+i)(21-i) = 421+i )
∴ 4tan-1(1/421) + 16tan-1(1/21) + 4tan-1(1/1985) + 3tan-1(1/239) = tan-1(1) (9-9)
ところが
tan-1(1/1985) + tan-1(1/421) = tan-1(3/1042) (9-10)
( ∵ (1985+i)(421+i) = 835684+i2406 = 802(1042+3i) )
∴ 16tan-1(1/21) + 3tan-1(1/239) + 4tan-1(3/1042) = π/4
（証明終了）
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/239) と tan-1(3/1042) を除いた部分を実行する．
(3) flg==2 で tan-1(1/21) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 247197[ミリ秒] = 4 分 07 秒 197 ミリ秒
(2) flg==1 のとき 117031[ミリ秒] = 1 分 57 秒 031 ミリ秒
(3) flg==2 のとき 131352[ミリ秒] = 2 分 11 秒 352 ミリ秒
これは一台の PC で実行した場合に驚くべきことにマチンの公式よりも約 38 秒ほども速く今まで述べたどの級数よりも速い．
この (9-1) の tan-1(3/1042) に正規化公式 (6-1) を適用して (a,b,±c）= (347,3,1) を代入すれば,
tan-1(3/1042) = tan-1(1/347) - tan-1(1/361577) (9-11)
となるので,これを (9-1) に代入すれば
◇ 16tan-1(1/21) + 3tan-1(1/239) + 4tan-1(1/347) - 4tan-1(1/361577) = π/4 (10-1)
が得られる．
これも以前のような方法で測定を試みた．
(0) 計算の要求精度を 10000（桁）とする．
(1) flg==0 で全てを実行する．
(2) flg==1 で tan-1(1/239) と tan-1(1/347) と tan-1(1/361577) を除いた部分を実行する．
(3) flg==2 で tan-1(1/21) を除いた部分を実行する．
結果は次のようであった．
(1) flg==0 のとき 270301[ミリ秒] = 4 分 30 秒 301 ミリ秒
(2) flg==1 のとき 117171[ミリ秒] = 1 分 57 秒 171 ミリ秒
(3) flg==2 のとき 153457[ミリ秒] = 2 分 33 秒 457 ミリ秒
これは一台の PC で実行した場合にマチンの公式よりも約 15 秒ほど速いが (9-1) より約 23 秒遅い．
(9-1) では分数の分子が 1 でないため級数の一項あたりの処理行数は増えるが tan-1 部は 3 部で良い，
一方 (10-1) では 分数の分子が 1 のため級数の一項あたりの処理行数は少ないが tan-1 部が一つ増えて 4 部になっている．
この微妙な駆引きによって結局 (10-1) の方が (9-1) より若干遅くなると考えられる．
このようなことが既存の公式の改良のヒントとなるかも知れない．
（ 運が良ければ (4-1) から (10-1) のどれかで私の新発見と言える公式があるかも知れない．）
さらにより本稿の趣旨に合致した次のような逆正接級数法の一般化が可能である．
§３.３.６ 単位複素弧度回転演算子法の計算結果を応用した逆正接関数のテーラー級数による円周率近似の一般的拡張
単位複素弧度回転演算子法の前提から，
◇ Mtan-1(S/C) + tan-1(C(M)/S(M)) = π/2 (1-1)
あるいは
◇ Mtan-1(1/H) + tan-1(C(M)/(1+S(M))) = π/4 (1-2)
が成り立つ．
あるいは別の言い方をすれば，適格な整数を U,V として，
◇ Mtan-1(1/H) = tan-1(V/U) ⇒ Mtan-1(1/H) + tan-1((U-V)/(U+V) = π/4 (1-3)
となる．
この公式から従来から良く知られていた公式を導いてみよう．
[1] H=1 のとき M=1 とすれば，
U=1,V=1,U+V=2,U-V=0 となるので，グレゴリーの公式
tan-1(1) = π/4 (2-1)
が導ける．
[2] H=2 のとき M=1 とすれば，
U=2,V=1,U+V=3,U-V=1 となるので，オイラーの公式
tan-1(1/2) + tan-1(1/3) = π/4 (2-2)
が導ける．
[3] H=3 のとき M=2 とすれば，
U=4,V=3,U+V=7,U-V=1 となるので，ハットンの公式
2tan-1(1/3) + tan-1(1/7) = π/4 (2-3)
が導ける．
[4] H=4 のとき M=3とすれば,
U=52,V=47,U+V=99,U-V=5 となるので，もう一つのハットンの公式
3tan-1(1/4) + tan-1(5/99) = π/4 (2-4)
が導ける．
[5] H=5 のとき M=4とすれば，
U=119,V=120,U+V=239,U-V=-1 となるので，マチンの公式
4tan-1(1/5) - tan-1(1/239) = π/4 (2-5-1)
が導ける．
（ これらは大変面白いことにここに掲げた限りにおいては後になるほど収束が早くなっている．
しかし，これ以降では PC で計算させる場合の処理時間の意味である限りにおいて，そうでもない．）
さらに (2-5-1) において，M=3とすれば,
U=55,V=37,(U+V,U-V)=2(46,9) となるので，
3tan-1(1/5) + tan-1(9/46) = π/4 (2-5-2)
tan-1(1/5) = tan-1(1/239) + tan-1(9/46) (2-5-3)
も成り立つのである．
ここで M=4, -tan-1(1/239) は先述の過越し商と過越し剰余にあたっている．
このように単位複素弧度回転演算子法はこれらの級数を導くための独立した入り口の役目を果たしているとも言えるのである．
この方法で侮れないのは，運が良ければマチンの公式より速い公式も得られることである．
例えば
H=28,M=22 とすれば
U=49324070120846121302260022403183,
V=49322325613363940973893167838056,
U+V=98646395734210062276153190241239,
U-V=1744507482180328366854565127,
(U+V)/(U-V))=56546.84588163495,
となる．これは U,V が極めて巨大である．
このときの私の PC での実行結果はは 10000 桁の円周率の計算において
全体の処理時間 220.990[秒]=3[分]40.990[秒]
主に tan-1(1/28) の処理時間 105.971[秒]=1[分]45.971[秒]
主に tan-1((U-V)/(U+V)) の処理時間 115.565[秒]=1[分]55.565[秒]
となった．これはマチンの公式より約 1[分]04[秒] 速い．
H=28,M=22 における 10000 桁の近似
これを更に改良して
22*ATN(1/28) + ATN(1/56546) + ATN((U-V)/(U+V)) = ATN(1) (2-5-4)
U-V =
-1475646841214443994950569897
U+V =
5578059094931149663647686662235665621
(U+V)/|U-V| =
3780077278.070448
処理時刻差 = 174486 ミリ秒
処理時間 = 2 分 54 秒 486 ミリ秒
が得られた．これはマチンの公式より約 1[分]50[秒] 速く,
私の PC 上では円周率の 10000 桁の計算で 3 分を切る唯一の方法である．(2010'09'23/Thu)
22*ATN(1/28)+ATN(1/56546)-ATN(|U-V|/(U+V))=ATN(1)による 10000 桁の近似
別の例として
H=14,M=11 とすれば
U=2947746030782,
V=2941761989627,
U+V=5889508020409,
U-V=5984041155,
(U+V)/(U-V))=984.2024591505337,
となる．
このときの私の PC での実行結果は 10000 桁の円周率の計算において
全体の処理時間 259.756[秒]=4[分]19.756[秒]
主に tan-1(1/14) に要する処理時間 134.098[秒] = 2[分]14.098[秒]
主に tan-1((U-V)/(U+V)) に要する処理時間 126.594[秒] = 2[分]06.594[秒]
となった．これもマチンの公式より約 25[秒] 速い．
これらの H,M の間では
M/H(22/28=11/14=0.785714285) ≒ π/4(=0.78539816…) (2-6)
が成り立っている．
級数法全体の傾向として
（１）1/H, (U-V)/(U+V) の値が小さいほど収束が速い．
（２）tan-1(x) の筋立てを複数本持つ場合は筋立てが少ない方が収束が速い．
それゆえ筋立てが 4 本もあるいわゆる高野の公式が最善の公式であるかどうかについて私は疑念を禁じえない．
案の定，高野の公式
◇ 12tan-1(1/49) + 32tan-1(1/57) - 5tan-1(1/239) + 12tan-1(1/110443) = π/4 (3)
によって私の PC で 10000 桁の円周率を計算してみたところ次のような結果となった．
高野の公式 による 10000 桁の近似
全体の処理時間 275.527[秒] = 4[分]35.527[秒]
主に tan-1(1/49) に要する処理時間 91.697[秒] = 1[分]31.697[秒]
主に tan-1(1/49) 以外に要する処理時間 183.441[秒] = 3[分]03.441[秒]
これはマチンの公式より約 10[秒] 速いにすぎない．
私も日本人だがマチンの公式より速い公式で自力で見つけた公式を 6 つは挙げることが出来る．
（ その内 5 つは高野の公式より速い．）
§３.４ 単位複素弧度回転演算子法の内で二冪商と三冪商との比較
単位複素弧度回転演算子法において，商 M が約 x 桁の１０進数であるとする．
このとき，さらに y, z を適当な正の整数として
M ≒ 10x ≒ 2y ≒ 3z (1-1)
であったとすれば，
y は二倍角公式のおよその実行回数であり，z は三倍角公式のおよその実行回数であると考えられる．
y ≒ x/log102 ≒ x/0.30103 ≒ 3.32x (1-2)
z ≒ x/log103 ≒ x/0.47712 ≒ 2.10x (1-3)
ここから一見すれば三倍冪の方が二倍冪より少ない演算回数で要求精度を満たすような期待を抱かせる．
ところが，二倍角公式を一回実行する場合の複素二項演算の回数は一回であるのに対し．
三倍角公式を一回実行する場合の複素二項演算の回数は倍角公式一回と加法定理一回の合計二回となる．
それで二冪商で実行した場合の複素二項演算の回数を N2,
三冪商で実行した場合の複素二項演算の回数を N3 とすれば，
N2 ≒ y ≒ x/log102 ≒ 3.32x (2-1)
N3 ≒ 2z ≒ 2x/log103 ≒ 2・2.10x ≒ 4.20x (2-2)
∴ N3/N2 ≒ 2log102/log103 ≒ 4.20/3.32 ≒ 1.265 (2-3)
結局，三冪商は二冪商と比較して約 1/4 ほど処理時間が長くなると予想される．
実際に 10000 桁の精度で円周率を求めた場合の私の例では，
三冪商で M = 3z = 310480 とした場合の処理時間は約 6979549 [msec]
二冪商で M = 2y = 216610 とした場合の処理時間は約 5255656 [msec]
その比を計算すれば，
6979549/5255656 ≒ 1.328
となって理論上の比 1.265 と 5% 程度の誤差の範囲で一致している．
（ 理論上の比とのズレが生じた原因は，
二冪商の計算プログラム内の二倍角公式の 16610 回の繰り返しループの命令文ブロックと比較して，
三冪商の計算プログラム内の三倍角公式の 10480 回の繰り返しループの命令文ブロックの主演算以外の補助処理が一行だけ多い
ことによると考えられる．）
同様に３を越える n 倍角公式の場合も回数を z' とすれば，
z' = x/log10n ＜ x/log103 ＜ x/log102
となるから一見少なくなったように見えるが，
一回あたりの n 倍角公式中に含まれる複素二項演算の回数が増える結果，全体の処理時間も増えるはずである．
ゆえに，単位複素弧度回転演算子法において，
勝手な n 倍角公式よりは二倍角公式を優先する方が処理時間が最短になると考えられる．
§３.７.１ 主としてオイラーの方法等に代表される逆正接関数の級数展開による円周率の近似計算公式に関して
判別式が２の二次拡大整数環の単数方程式を応用した逆正接公式の一般的生成法
この件に関して基礎的な事柄は私のホーム頁の冒頭の目次から
（１）「 判別式が２の二次拡大整数環の基本性質 」
（２）「 整数の三平方の定理の直交する二辺の差 」
の節で述べた内容とも一部が重複せざるを得ない．
それゆえ上記の節を再読する手間を省く為に基本的な定義の部分を簡単に繰り返しておくこととする．
先ず，ここで私が判別式が２の二次拡大単数方程式と呼んでいるのは次のものである．
n を順序数, fn,gn を適当な二つ一組の正の有理整数として，
fn + gn√2 = ( 1 + √2 )n ( n = 1,2,3,… ) (1-1-1)
と定義する．
するとこれらは n の奇数 ( n ≡ 1 (mod 2)) か偶数 ( n ≡ 0 (mod 2)) に応じて，
fn2 - 2gn2 = -1 ( n ≡ 1 (mod 2)) (1-1-2)
fn2 - 2gn2 = 1 ( n ≡ 0 (mod 2)) (1-1-3)
に分かれる．
さらに三つ一組の有理整数を xn,yn,zn として，
xn = ( fn + 1 )/2, yn = ( fn - 1 )/2, zn = gn (1-2-1)
と定義する．
するとこれらの xn,yn,zn の間には，
zn ≡ n (mod 2) (1-2-2)
xn - yn = 1 (1-2-3)
xn2 + yn2 = zn2 ( n ≡ 1 (mod 2)) (1-2-4-1)
xn2 + yn2 = zn2 + 1 ( n ≡ 0 (mod 2)) (1-2-4-2)
などの関係がある．
ここで (1-2-3-1) はお馴染みの三平方の定理（ピタゴラスの定理）に属している．
しかし (1-2-3-2) は一余り三平方の定理とでも呼ぶべきものとなっている．
これらの序列の初めの方の部分を数表として掲げてみると，
[数表１] 判別式が２の二次拡大単数方程式から求められた fn,gn,xn,yn,zn の値（√2単数表）
n fn gn xn yn zn
1 1 1 1 0 1
2 3 2 2 1 2
3 7 5 4 3 5
4 17 12 9 8 12
5 41 29 21 20 29
6 99 70 50 49 70
7 239 169 120 119 169
8 577 408 289 288 408
9 1393 985 697 696 985
10 3363 2378 1682 1681 2378
これらの数表からも伺えることだが（この数表のことを以後 √2単数表と呼ぶことにしよう．）
fn+1 = fn + 2gn, gn+1 = fn + gn (1-3-1-1)
fn+2 = 2fn+1 + fn, gn+2 = 2gn+1 + gn (1-3-1-2)
などが成り立つ．
さて更に，√2 の近似に関して，
fn/gn ＜ √2 ＜ 2gn/fn ( n ≡ 1 (mod 2 )) (1-3-2-1)
2gn/fn ＜ √2 ＜ fn/gn ( n ≡ 0 (mod 2 )) (1-3-2-2)
が成り立つ．
ゆえにまた，
1 + fn/gn ＜ 1 + √2 ＜ 1 + 2gn/fn ( n ≡ 1 (mod 2 )) (1-3-3-1)
1 + 2gn/fn ＜ 1 + √2 ＜ 1 + fn/gn ( n ≡ 0 (mod 2 )) (1-3-3-2)
ところが，(1-3-1-1),(1-3-1-2) および
tan(π/8) = 1/( 1 + √2 ) (1-3-4)
から，
fn/fn+1 ＜ tan(π/8) ＜ gn/gn+1 ( n ≡ 1 (mod 2 )) (1-3-5-1)
gn/gn+1 ＜ tan(π/8) ＜ fn/fn+1 ( n ≡ 0 (mod 2 )) (1-3-5-2)
それゆえ，また，
tan-1(fn/fn+1) ＜ π/8 ＜ tan-1(gn/gn+1) ( n ≡ 1 (mod 2 )) (1-3-6-1)
tan-1(gn/gn+1) ＜ π/8 ＜ tan-1(fn/fn+1) ( n ≡ 0 (mod 2 )) (1-3-6-2)
これはこの公式だけでも十分面白いが，しかし，まだ不等式であるにすぎない．
実はこれを等式に変える工夫も存在する．すなわち次の定理が成り立ち，証明も可能である．
[定理]
tan-1(fn/fn+2m-1) + tan-1(gn/gn+2m-1) = tan-1(1/f2m-1) (1-4-1)
tan-1(fn/fn+2m-1) - tan-1(gn/gn+2m-1) = tan-1(1/(f2m-1+(-1)n4gngn+2m-1)) (1-4-2)
[証明]
証明を簡単にするために，新たに F,G,f,g,f*,g* を用いて
F = fn+2m-1, G = gn+2m-1, f = fn, g = gn, f* = f2m-1, g* = g2m-1 であるとする．
定理の前提から，
F + G√2 = ( f + g√2 )( f* + g*√2 )
∴ F = ff* + 2gg*, G = fg* + f*g (1-5-1)
ここで, 新たに U,V を用いて，次のような複素数の乗法を考える
U + iV = ( F + if )( G + ig )
すると，
U = FG - fg, V = Fg + fG (1-5-2)
∴ U = ( ff* + 2gg* )( fg* + f*g ) - fg
= f2f*g* + ff*2g + 2fgg*2 + 2f*g2g* - fg
= f*g*( f2 + 2g2 ) + fg( f*2 + 2g*2 - 1 )
= f*g*( 2f2 - (-1)n ) + 2fgf*2
( ∵ f2 - 2g2 = (-1), f*2 - 2g*2 = -1 )
= -(-1)nf*g* + 2ff*( fg* + f*g )
= -(-1)nf*g* + 2ff*G
= f*( -(-1)ng* + 2fG ) (1-5-3)
V = ( ff* + 2gg* )g + f( fg* + f*g )
= 2ff*g + 2g2g* + f2g*
= 2f( f*g + fg* ) - ( f2 - 2g2 )g*
= -(-1)ng* + 2fG (1-5-4)
∴ V/U = 1/f* (1-5-5)
∴ ( F + if )( G + ig )/(-(-1)ng*+2fG) = f* + i (1-5-6)
これで (1-4-1) が証明された．
再び U,V を，
U + iV = ( F + if )( G - ig ) (1-6-1)
とすれば，
U = FG + fg, V = fG - Fg (1-6-2)
U = ( ff* + 2gg* )( fg* + f*g )
= f2f*g* + ff*2g + 2fgg*2 + 2f*g2g* + fg
= f*g*( f2 + 2g2 ) + fg( f*2 + 2g*2 + 1 )
= f*g*( 4g2 +(-1)n ) + 4fgg*2
( ∵ f2 - 2g2 = (-1), f*2 - 2g*2 = -1 )
= (-1)nf*g* + 4gg*( f*g + fg* )
∴ U = (-1)nf*g* + 4gg*G (1-6-3)
V = f( fg* + f*g ) - ( ff* + 2gg* )g
= f2g* + ff*g - ff*g - 2g2g*
= g*( f2 - 2g2 )
= (-1)ng* (1-6-4)
∴ V/U = 1/( f* + (-1)n4gG ) (1-6-5)
∴ ( F + if )( G - ig )/((-1)ng*) = ( f* + (-1)n4gG ) + i (1-6-6)
これで (1-4-2) が証明された．
（証明終了）
さらに (1-4-1),(1-4-2) から，
[定理]
2tan-1(fn/fn+2m-1) - tan-1(1/(f2m-1+(-1)n4gngn+2m-1)) = tan-1(1/f2m-1) (1-4-3)
2tan-1(gn/gn+2m-1) + tan-1(1/(f2m-1+(-1)n4gngn+2m-1)) = tan-1(1/f2m-1) (1-4-4)
が成り立つ．
さらに n = 1 のときは，
f1 = g1 = 1
であることに注意すれば，
[定理]
tan-1(1/f2m) + tan-1(1/g2m) = tan-1(1/f2m-1) (1-7-1)
tan-1(1/f2m) - tan-1(1/g2m) = tan-1(1/(f2m-1-4g2m)) = - tan-1(1/f2m+1) (1-7-2)
2tan-1(1/f2m) = tan-1(1/f2m-1) - tan-1(1/f2m+1) (1-7-3)
2tan-1(1/g2m) = tan-1(1/f2m-1) + tan-1(1/f2m+1) (1-7-4)
となる．
[検討]
n = 1 のとき (1-4-2) の右辺の分母は次のように変形できる．
f2m-1 + (-1)n4gngn+2m-1
= f2m-1 - 4g1g2m
= f2m-1 - 4g2m
= f2m-1 - 2(f2m-1+g2m-1) - 2g2m
= -f2m-1 - 2g2m-1 - 2g2m
= -f2m-1 - 2(g2m-1+g2m)
= -f2m-1 - 2f2m
= -f2m+1
( ∵ f2m+1 = f2m-1 + 2f2m ) //
さらに m = 1 のときは
[定理]
tan-1(fn/fn+1) + tan-1(gn/gn+1) = tan-1(1) (1-8-1)
tan-1(fn/fn+1) - tan-1(gn/gn+1) = (-1)ntan-1(1/f2n+1) (1-8-2)
2tan-1(fn/fn+1) - (-1)ntan-1(1/f2n+1) = tan-1(1) (1-8-3)
2tan-1(gn/gn+1) + (-1)ntan-1(1/f2n+1) = tan-1(1) (1-8-4)
[検討]
m = 1 のとき (1-4-2) の右辺の分母は次のように変形できる．
f2m-1 + (-1)n4gngn+2m-1
= f1 + (-1)ngn+1gn
= 1 + (-1)n4gngn+1
= (-1)n( (-1)n + 4gngn+1 )
= (-1)n( fn2 - 2gn2 + 4gngn+1 )
= (-1)n( (fn+1-2gn)fn - (fn+1-fn)gn + 4gn+1gn )
= (-1)n( fn+1fn - fngn - fn+1gn + 4gn+1gn )
= (-1)n( fn+1fn - (fn+1+fn)gn + 4gn+1gn )
= (-1)n( fn+1fn - 2gn+1gn + 4gn+1gn )
= (-1)n( fn+1fn + 2gn+1gn )
= (-1)nf2n+1
( ∵ f2n+1 = fn+1fn + 2gn+1gn ) //
これらの無限通りの公式の生成法は私の独自の発見であり研究の結果である．
これらの一連の公式たちから少なくとも言える事は
φ(x) + φ(y) = φ(1) = π/4 ( x,y は有理数 ) (1-9)
を満たすような適当な写像 φ の一つとして
tan-1 を挙げることができ，
無限に多くの有理数解 (x,y) を含むということである．
通常の水平な実数直線ではこれらのことは適わない．
この様な意味において tan-1 は正に "丸い地球の水平線" に当たっている．
私は，これら一連の方法を飛燕奥義八分円逆正接法と命名した．（ 略称 飛燕矢車 ）
ここでも 飛燕六三四の場合と同様に √2 単数表すなわち判別式が２の双曲線が使えるのである．
§３.７.１ √2単数表を応用した逆正接等式の具体例
前節の公式 (1-4-1) 等の具体的な応用を調べてみよう.
ただしここでは級数の収束の遅速は問題にせず，ただ数学的に不正な公式でなければそれで良いとしよう．
m = 1 を前提として，n = 1 から順次 n を増加させると次のような公式が得られる.
（１）n = 1 のとき (1-8-1) において
f1 = 1, f2 = 3, g1 = 1, g2 = 2 を代入すれば,
tan-1(1/3) + tan-1(1/2) = tan-1(1) = π/4 (1-1-1)
これは正にオイラーの公式以外の何者でもない．それゆえ前節で述べた方法はオイラーの公式の自然な拡張である．
さらに (1-8-3) で f2n+1 = f3 = 7 を代入すれば,
2tan-1(1/3) + tan-1(1/7) = π/4 (1-1-2)
を得る．これはハットンの第一公式である．
さらに (1-8-4) から,
2tan-1(1/2) - tan-1(1/7) = π/4 (1-1-3)
が得られる．この公式はハットンによって発見されていたようである．
以下はこのような要領に従うものとして初めの数個について結果だけを列挙しておこう．
（２）n = 2 のとき (1-8-1) において
f2 = 3, f3 = 7, g2 = 2, g3 = 5 を代入すれば,
tan-1(3/7) + tan-1(2/5) = π/4 (1-2-1)
さらに (1-8-3),(1-8-4) で f2n+1 = f5 = 41 を代入すれば,
2tan-1(3/7) - tan-1(1/41) = π/4 (1-2-2)
2tan-1(2/5) + tan-1(1/41) = π/4 (1-2-3)
（３）n = 3 のとき (1-8-1) において
f3 = 7, f4 = 17, g3 = 5, g4 = 12 を代入すれば,
tan-1(7/17) + tan-1(5/12) = π/4 (1-3-1)
さらに (1-8-3),(1-8-4) で f2n+1 = f7 = 239 を代入すれば,
2tan-1(7/17) + tan-1(1/239) = π/4 (1-3-2)
2tan-1(5/12) - tan-1(1/239) = π/4 (1-3-3)
このときさらに (1-3-3) では
tan-1(5/12) = 2tan-1(1/5) (1-3-4)
であることに注意すれば,
4tan-1(1/5) - tan-1(1/239) = π/4 (1-3-5)
を得る．これはマチンの公式である．
ここはまた前掲の √2単数表の中では稀な箇所である，というのは,
g3 = 5, g4 = 12
であるのだが，この箇所だけが
√(g32+g42) = √(52+122) = 13 (1-3-6)
と成って右辺の平方根が開けるのである．
これ以外で,
gn+12 + gn2 = k2, ( k は整数 ) (1-3-7)
が成り立つところはこの数表の無限に及ぶ全体の中でここだけであろうと予想されるが残念ながら私には証明出来ない．
とにかくマチンの公式はその稀な箇所にある不思議な公式なのである．（ そして実用的にも簡単な公式の割りに妙に収束が速い．）
（４）n = 4 のとき (1-8-1) において
f4 = 17, f5 = 41, g3 = 12, g4 = 29 を代入すれば,
tan-1(17/41) + tan-1(12/29) = π/4 (1-4-1)
さらに (1-8-3),(1-8-4) で f2n+1 = f9 = 1393 を代入すれば,
2tan-1(17/41) - tan-1(1/1393) = π/4 (1-4-2)
2tan-1(12/29) + tan-1(1/1393) = π/4 (1-4-3)
（５）n = 5 のとき (1-8-1) において
f5 = 41, f6 = 99, g5 = 29, g6=70 を代入すれば,
tan-1(41/99) + tan-1(29/70) = π/4 (1-5-1)
さらに (1-8-3),(1-8-4) で f2n+1 = f11 = 8119 を代入すれば,
2tan-1(41/99) + tan-1(1/8119) = π/4 (1-5-2)
2tan-1(29/70) - tan-1(1/8119) = π/4 (1-5-3)
§３.７.２ √2 単数表を応用した類似の逆正接関数の等式
ここでは前節と類似の逆正接関数の等式の内でこの数表の特徴的な公式を幾つか挙げてみよう．
便宜のためここに数表を再掲しておく．
[数表１] √2 単数表（再掲）
n fn gn xn yn zn
1 1 1 1 0 1
2 3 2 2 1 2
3 7 5 4 3 5
4 17 12 9 8 12
5 41 29 21 20 29
6 99 70 50 49 70
7 239 169 120 119 169
8 577 408 289 288 408
9 1393 985 697 696 985
10 3363 2378 1682 1681 2378
次の公式が成り立つ．
[公式]
tan-1(gn/fn) + tan-1(gn-1/gn+1) = tan-1(1) (n≧2) (1-1)
tan-1(fn/2gn) + tan-1(fn-1/fn+1) = tan-1(1) (n≧2) (1-2)
[証明]
gn-1 = fn-gn, gn+1 = fn+gn
となるので（以後の証明では記号を簡単にするため fn = f,gn = g とする．）
(fn+ign)(gn+1+ign-1)
= (f+ig)((f+g)+i(f-g)) = (f2+fg-fg+g2)+i(f2-fg+fg+g2)
= (f2+g2)(1+i)
ゆえに (1-1) が示された．
fn-1 = 2gn-fn, fn+1 = fn+2gn
となるので
(2gn+ifn)(fn+1+ifn-1)
= (2g+if)((f+2g)+i(2g-f))
= (2fg+4g2-2fg+f2)+i(4g2-2fg+f2+2fg)
= (f2+4g2)(1+i)
ゆえに (1-2) が示された．
（証明終了）
(1-1) の例を列挙すれば，
(1) n=2
tan-1(2/3) + tan-1(1/5) = tan-1(1)
(2) n=3
tan-1(5/7) + tan-1(1/6) = tan-1(1)
(3) n=4
tan-1(12/17) + tan-1(5/29) = tan-1(1)
(4) n=5
tan-1(29/41) + tan-1(6/35) = tan-1(1)
(5) n=6
tan-1(70/99) + tan-1(29/169) = tan-1(1)
(1-2) の例を列挙すれば，
(1) n=2
tan-1(3/4) + tan-1(1/7) = tan-1(1)
(2) n=3
tan-1(7/10) + tan-1(7/17) = tan-1(1)
(3) n=4
tan-1(17/24) + tan-1(7/41) = tan-1(1)
(4) n=5
tan-1(41/58) + tan-1(17/99) = tan-1(1)
(5) n=6
tan-1(99/140) + tan-1(41/239) = tan-1(1)
さらに √2 単数表の特性として次の公式が成り立つことも面白い．
[公式]
2tan-1(fn/fn+1) = tan-1(y2n+1/x2n+1) ( n≡1 (mod 2)) (2-1-1)
2tan-1(fn/fn+1) = tan-1(x2n+1/y2n+1) ( n≡0 (mod 2)) (2-1-2)
2tan-1(gn/gn+1) = tan-1(x2n+1/y2n+1) ( n≡1 (mod 2)) (2-2-1)
2tan-1(gn/gn+1) = tan-1(y2n+1/x2n+1) ( n≡0 (mod 2)) (2-2-2)
[証明]
(fn+1+ifn)2 = (fn+1)2-(fn)2+i2fn+1fn
= (fn+2gn+fn)(fn+2gn-fn)+i2(fn+2gn)fn
= 2(fn+gn)gn+i(2fn2+2fngn)
= (2fngn+gn2)+i(2fn2+2fngn)
= (fn2+4fngn+2gn2-(fn2-2gn2))/2+i(fn2+4fngn+2gn2+(fn2-2gn2))/2
= (f2n+1-(-1)n)+i(f2n+1+(-1)n)
(∵ fn2+4fngn+2gn2=f2n+1,fn2-2gn2=(-1)n)
= x2n+1+iy2n+1 ( n≡1 (mod 2))
または
= y2n+1+ix2n+1 ( n≡0 (mod 2))
(∵ x2n+1=(f2n+1+1)/2,y2n+1=(f2n+1-1)/2
ゆえに (2-1-1),(2-1-2) が示された．
(gn+1+ign)2 = (gn+1)2-(gn)2+i2gn+1gn
= (fn+gn+gn)(fn+gn-gn)+i2(fn+gn)gn
= (fn+2gn)fn+i(2gn2+2fngn)
= (2fngn+fn2)+i(2gn2+2fngn)
= (fn2+4fngn+2gn2+(fn2-2gn2))/2+i(fn2+4fngn+2gn2-(fn2-2gn2))/2
= (f2n+1+(-1)n)+i(f2n+1-(-1)n)
= y2n+1+ix2n+1 ( n≡1 (mod 2))
または
= x2n+1+iy2n+1 ( n≡0 (mod 2))
ゆえに (2-2-1),(2-2-2) が示された．
（証明終了）
これらについても序列の初めの方の部分だけを列挙しておこう．
(1) n=1
2tan-1(1/3) = tan-1(3/4)
2tan-1(1/2) = tan-1(4/3)
(2) n=2
2tan-1(3/7) = tan-1(21/20)
2tan-1(2/5) = tan-1(20/21)
(3) n=3
2tan-1(7/17) = tan-1(119/120)
2tan-1(5/12) = tan-1(120/119)
(4) n=4
2tan-1(17/41) = tan-1(697/696)
2tan-1(12/29) = tan-1(696/697)
(5) n=5
2tan-1(41/99) = tan-1(4059/4060)
2tan-1(29/70) = tan-1(4060/4059)
これらの公式の右辺は √2 単数表の奇数行だけが対応している．
そこで偶数行が対応するような類似の公式として，次のような公式が挙げられる．
[公式]
（ √2 単数表に n=0,f0=1,g0=0 を予め追加しておくものとする．）
tan-1(fn-1/fn)+tan-1(fn/fn+1) = tan-1(x2n/y2n) ( n≡1 (mod 2)) (3-1-1)
tan-1(fn-1/fn)+tan-1(fn/fn+1) = tan-1(y2n/x2n) ( n≡0 (mod 2)) (3-1-2)
tan-1(gn-1/gn)+tan-1(gn/gn+1) = tan-1(y2n/x2n) ( n≡1 (mod 2)) (3-2-1)
tan-1(gn-1/gn)+tan-1(gn/gn+1) = tan-1(x2n/y2n) ( n≡0 (mod 2)) (3-2-2)
[証明]
(fn+ifn-1)(fn+1+ifn)
= (fn+i(2gn-fn))((fn+2gn)+ifn)
= fn2+2fngn-2fngn+fn2+i(fn2+4gn2-fn2)
= 2fn2+i4gn2
= (f2n+(-1)n)+i(f2n-(-1)n)
(∵ f2n=fn2+2gn2,f2n2-2g2n2=(-1)n)
= tan-1(x2n/y2n) ( n≡1 (mod 2))
または
= tan-1(y2n/x2n) ( n≡0 (mod 2))
(∵ x2n=(f2n+1)/2,y2n=(f2n-1)/2)
これで (3-1-1),(3-1-2) が示された．
(gn+ign-1)(gn+1+ign)
= (gn+i(fn-gn))((fn+gn)+ign)
= fngn+gn2-fngn+gn2+i(gn2+fn2-gn2)
= 2gn2+ifn2
= (f2n-(-1)n)/2+i(f2n+(-1)n)/2
= tan-1(y2n/x2n) ( n≡1 (mod 2))
または
= tan-1(x2n/y2n) ( n≡0 (mod 2))
これで (3-2-1),(3-2-2) が示された．
（証明終了）
これらについても序列の初めの方の部分だけを列挙しておこう．
(1) n=1
tan-1(1/1)+tan-1(1/3) = tan-1(2/1)
tan-1(0/1)+tan-1(1/2) = tan-1(1/2)
(2) n=2
tan-1(1/3)+tan-1(3/7) = tan-1(8/9)
tan-1(1/2)+tan-1(2/5) = tan-1(9/8)
(3) n=3
tan-1(3/7)+tan-1(7/17) = tan-1(50/49)
tan-1(2/5)+tan-1(5/12) = tan-1(49/50)
(4) n=4
tan-1(7/17)+tan-1(17/41) = tan-1(288/289)
tan-1(5/12)+tan-1(12/29) = tan-1(289/288)
(5) n=5
tan-1(17/41)+tan-1(41/99) = tan-1(1682/1681)
tan-1(12/29)+tan-1(29/70) = tan-1(1681/1682)
また次のような公式も成立する．
[公式]
tan-1(fn-1/fn) - tan-1(fn/fn-1) = -(-1)ntan-1(1/g2n) (4-1-1)
tan-1(gn-1/gn) - tan-1(gn/gn-1) = (-1)ntan-1(1/g2n) (4-1-2)
（証明のヒント：g2n=2fngn,fn2-2gn2=(-1)n ）
tan-1(fn-1/fn) - tan-1(fn/fn-1) + (-1)ntan-1(1/f2n) = tan-1(1) (4-2-1)
tan-1(gn-1/gn) - tan-1(gn/gn-1) - (-1)ntan-1(1/f2n) = tan-1(1) (4-2-2)
2tan-1(fn-1/fn) + (-1)ntan-1(1/f2n) + (-1)ntan-1(1/g2n) = tan-1(1) (4-3-1)
2tan-1(gn-1/gn) - (-1)ntan-1(1/f2n) - (-1)ntan-1(1/g2n) = tan-1(1) (4-3-2)
これらの例の若干を挙げておけば
tan-1(1/17) + tan-1(1/12) = tan-1(1/7)
tan-1(1/99) + tan-1(1/70) = tan-1(1/41)
tan-1(1/577) + tan-1(1/408) = tan-1(1/239)
2tan-1(1/3) + tan-1(1/17) + tan-1(1/12) = tan-1(1)
2tan-1(3/7) - tan-1(1/99) - tan-1(1/70) = tan-1(1)
2tan-1(2/5) + tan-1(1/99) + tan-1(1/70) = tan-1(1)
などがある．
§ ω. 終わりに代えて．
1949 年に世界初の汎用コンピュータ ENIAC が 70 時間かけて π を 2037 桁計算した．
真空管式でプログラム内臓式ではなく，重量は 30 トンもあり消費電力は 150 kW で
倉庫一つ分ほどの大きさもあったそうである．
それから約 62 年ほど経って今では
半導体メモリでプログラム内蔵式で，重量は 2.5 Kg 未満，消費電力 20 w 未満で
膝の上に乗る程度の PC であっても
（ 実は私はこの PC をプロバイダ様の乗り換え時の条件で僅か 500 円で入手したのだが ^^ ）
そんな誰でもが最低限度持っていそうな PC であっても
わずか 5 分にも満たない間に 10000 桁ほどの π を計算してしまうのである．
時代の進歩とは実にも恐ろしきものである．
円周率は適えども 己が命数は 数え難し (活己)
( 2011'09'04_Sun )
本稿終わり．
