- πに近い分数として知られる22/7や355/113は偶然の産物ではなく、連分数展開から導かれる収束分数です。
- 連分数とは「整数部分を取り出し、残りを逆数にする」操作を繰り返して実数を表す方法で、途中で止めると有理近似が得られます。
- 355/113が突出して精度が高い理由は、連分数の次の項が292という大きな数であるため、次の記録更新まで分母が33102まで跳び上がるからです。
1. 「πに近い分数」とは？
22/7は、3.1428… と円周率に近い分数です。
さらに近い分数としては、5世紀の中国の数学者 祖冲之が 355/113 を見つけています1。
pi
;=> 3.141592653589793d0
(coerce (/ 22 7) 'double-float)
;=> 3.142857142857143d0
(coerce (/ 355 113) 'double-float)
;=> 3.1415929203539825d0Code language: Lisp (lisp)22/7は、小数第2位まで、335/113は、小数第6位まで同じ数字になっています。
これらの分数は、「たまたま」近い値なのでしょうか。
実は、この2つには共通の数学的な背景があります。
πを有理数で近似するという問いをCommon Lispで掘り下げてみます。
1.1. 円周率 π に近い有理数は無限にある
有理数とは、分子と分母がともに整数である数のことです。
πは無理数で、有理数では表せません2。
しかし、πに近い有理数は、いくらでも作れます3。
たとえば、πの小数展開を途中で切るだけで、どんどん近づいていくことがわかります。
3141592653589793 / 1000000000000000分母をもっと大きくすれば、さらに近い分数が作れます。
「最も近いもの」という唯一の答えはないのです。
1.2. 記録更新列という考え方
そこで考えるのが、「分母が小さいわりに、特によく近似している分数はどれか」という問題です。
πに近い分数を、分母を1から順に増やしていき、そのつど「今の時点での最良近似」を更新していく、という考え方で探してみましょう。
(defun record-breakers-data (x max-den)
"x の最良近似のうち、記録を更新した (p q err) のリストを返す。"
(let ((best-err most-positive-double-float)
(result '()))
(loop for q from 1 to max-den do
(let* ((p (round (* x q)))
(err (abs (- x (/ p q)))))
(when (< err best-err)
(setf best-err err)
(push (list p q err) result))))
(nreverse result)))Code language: Lisp (lisp)(let* ((p (round (* x q))))) では、x に対して各分母 q に対する最も近い分子 p を求めています4。
これを元に、誤差が更新されたときだけ記録します。
次のように実行すると、記録を塗り替えた分数が順に出てきます。
(let ((data (record-breakers-data pi 200000)))
(dolist (row data)
(destructuring-bind (p q err) row
(format t "~A/~A error=~,4E~%" p q err))))Code language: Lisp (lisp)3/1 error=1.4159E-1
22/7 error=1.2645E-3
333/106 error=9.4175E-5
355/113 error=2.6676E-7
103993/33102 error=5.7789E-10
104348/33215 error=3.3154E-10
...数字を眺めると、355/113で誤差がぐっと小さくなったことがわかります。
1.3. 両対数グラフで誤差の減り方を見る
両対数グラフで誤差の減り方を描いてみましょう。
まず vgplot をロードします。
vgplotはgnuplotをバックエンドにしてCommon Lispからグラフを描くライブラリです5。
(ql:quickload :vgplot)Code language: Lisp (lisp)あとは、プロットしてみます。
(defun plot-error-loglog (max-den)
(let* ((data (record-breakers-data pi max-den))
(xs (mapcar #'second data))
(ys (mapcar #'third data)))
(vgplot:loglog xs ys ";error;with linespoints pt 7")
(vgplot:title "Error in best rational approximations to pi")
(vgplot:xlabel "denominator q")
(vgplot:ylabel "|pi - p/q|")))
(plot-error-loglog 10000000)Code language: Lisp (lisp)グラフ全体はおおむね直線になりますが、355/113の点だけ大きく下に外れます。
分母113のわりに誤差が格段に小さく、次の記録更新（分母33102）との間に大きな段差があることが視覚的に分かります6。
まとめると、記録を塗り替えた分数はこうなります。
| 分数 | 値 | πとの誤差 |
|---|---|---|
| 3/1 | 3.0 | 0.14159… |
| 22/7 | 3.14285… | 0.00126… |
| 333/106 | 3.14150… | 0.0000942… |
| 355/113 | 3.14159292… | 0.000000266… |
| 103993/33102 | 3.14159265301… | 5.78×10⁻¹⁰ |
2. 連分数が答えを持っている
実は、これらの分数の求め方には、もう一つの方法があります。
それが、「連分数」です。
連分数（continued fraction）とは、数を次の形で表す方法です7。
x = a₀ + 1 / (a₁ + 1 / (a₂ + 1 / (a₃ + ...)))「整数部分を取り出し、残りを逆数にする」を繰り返す、というアルゴリズムで求めることができます。
連分数は、ある実数 x に対して、+ 1/a で調整しながら近似していく方法です。
(defun continued-fraction-terms (x n)
"x の連分数の最初の n 項を返す。"
(loop with y = (coerce x 'double-float)
repeat n
collect (let ((a (floor y)))
(setf y (- y a))
(if (zerop y)
a
(prog1 a (setf y (/ 1.0d0 y)))))))
(defun convergents-from-cf (terms)
"連分数の項列から収束分数 (p q) のリストを返す。"
(let ((p-2 0) (p-1 1)
(q-2 1) (q-1 0)
(result '()))
(dolist (a terms (nreverse result))
(let ((p (+ (* a p-1) p-2))
(q (+ (* a q-1) q-2)))
(push (list p q) result)
(setf p-2 p-1 p-1 p
q-2 q-1 q-1 q)))))Code language: Lisp (lisp)πを連分数で展開するとこう書けます。
x₀ = π ≈ 3.14159...
a₀ = floor(x₀) = 3
x₁ = 1 / (x₀ - 3) ≈ 7.0625...
a₁ = floor(x₁) = 7
x₂ = 1 / (x₁ - 7) ≈ 15.9966...
a₂ = floor(x₂) = 15
x₃ = 1 / (x₂ - 15) ≈ 1.0034...
a₃ = floor(x₃) = 1
x₄ = 1 / (x₃ - 1) ≈ 292.6345...
a₄ = floor(x₄) = 292 ← ここが大きいつまり、πの連分数は、どこまでも続きます。
π = 3 + 1 / (7 + 1 / (15 + 1 / (1 + 1 / (292 + ...))))連分数には省略記法があり、以下のように書きます。
π = [3; 7, 15, 1, 292, 1, 1, 1, 2, ...]2.1. 収束分数を出す
この連分数はどこまでも続きますが、途中で止めれば有理数になります。
[3] = 3/1
[3; 7] = 3 + 1/7 = 22/7
[3; 7, 15] = 3 + 1/(7 + 1/15) = 333/106
[3; 7, 15, 1] = 3 + 1 / (7 + 1 / (15 + 1 / (1))) = 355/113これを収束分数（convergent）といいます。
記録更新列の主要な項は、この収束分数列と一致しています。
連分数の係数（各段の整数部分）と、収束分数を求めます。
収束分数の漸化式は p_n = a_n * p_{n-1} + p_{n-2} という形で、前の2項から次を計算します。
初期値は p_{-2}=0, p_{-1}=1, q_{-2}=1, q_{-1}=0 です。
(let* ((terms (continued-fraction-terms pi 6))
(convs (convergents-from-cf terms)))
(format t "terms: ~A~%" terms)
(dolist (pq convs)
(destructuring-bind (p q) pq
(format t "~A/~A~%" p q))))Code language: Lisp (lisp)terms: (3 7 15 1 292 1)
3/1
22/7
333/106
355/113
103993/33102
104348/33215つまり、355/113がよい近似である理由は、収束分数に含まれているから、と言えます。
また、この連分数の次の項は 292 という大きな数です8。
そのため、355/113の次に精度が上がるまでに、分母は33102になってしまいます。
2.2. 近似値の収束をグラフで見る
「分母 q と近似値 p/q」のグラフを描きます。
(defun plot-approximations (max-den)
(let* ((data (record-breakers-data pi max-den))
(xs (mapcar #'second data))
(ys (mapcar (lambda (row)
(destructuring-bind (p q err) row
(declare (ignore err))
(/ p q)))
data))
(pi-ys (make-list (length xs) :initial-element pi)))
(vgplot:plot xs ys ";p/q;with linespoints pt 7"
xs pi-ys ";pi;with lines")
(vgplot:title "Rational approximations approaching pi")
(vgplot:xlabel "denominator q")
(vgplot:ylabel "p/q")))
(plot-approximations 500)Code language: Lisp (lisp)p/qの点が上下を行き来しながらπに収束していきます。
収束分数が交互に大きい側と小さい側を取ることが視覚的に分かります。
3. 誤差の収束
この連分数を使うと、最良近似を効率的に求めることができます。
(defun record-breakers-data/convergents (x max-den)
"x の最良近似のうち、記録を更新した (p q err) のリストを、連分数の収束分数と中間収束分数を使って返す。"
(labels ((emit (p q best-err result)
(let ((err (abs (- x (/ p q)))))
(if (< err best-err)
(values err (cons (list p q err) result))
(values best-err result)))))
(let* ((y (coerce x 'double-float))
(p-2 0) (p-1 1)
(q-2 1) (q-1 0)
(best-err most-positive-double-float)
(result '()))
(loop
(let* ((a (floor y))
(p (+ (* a p-1) p-2))
(q (+ (* a q-1) q-2)))
(when (> q max-den)
(return (nreverse result)))
(when (> a 1)
(loop for k from 1 below a do
(let ((pp (+ p-2 (* k p-1)))
(qq (+ q-2 (* k q-1))))
(when (> qq max-den)
(return))
(multiple-value-setq (best-err result)
(emit pp qq best-err result)))))
(multiple-value-setq (best-err result)
(emit p q best-err result))
(let ((frac (- y a)))
(when (zerop frac)
(return (nreverse result)))
(setf y (/ 1.0d0 frac)
p-2 p-1 p-1 p
q-2 q-1 q-1 q)))))))Code language: Lisp (lisp)ここで、両対数グラフで誤差をもう一度見てみます。
このグラフは、直線 y = C / q^2 に概ね重なります。
これにも、理論的な裏付けがあります。
ディリクレの近似定理では、良い近似における誤差は |π - p/q| < 1/q² を満たすことが知られており9、両辺の対数を取るとグラフの傾きはだいたい -2 になります。
ただし、π では連分数にしばしば 292 のような突出した項があるので、355/113 付近だけ誤差が異様に小さくなります。
それがグラフの飛び出した点として現れます。
3.1. まとめ
連分数の性質を実際にコードで動かして図にすると、「355/113だけ突出して良い」という事実が直感として入ってきます。
これが計算で探索することの面白さだと思います。
Common LispはREPLで少しずつ試しながら関数を組み合わせ、図を出すという作業に向いています。
その手触りが、数学の構造を感覚として掴む助けになります。
;; 最小構成で試すなら
(ql:quickload :vgplot)
(defun plot-pi (max-den)
(let* ((best-err most-positive-double-float)
(xs '()) (ys '()))
(loop for q from 1 to max-den do
(let* ((p (round (* pi q)))
(err (abs (- pi (/ p q)))))
(when (< err best-err)
(setf best-err err)
(push q xs)
(push err ys))))
(vgplot:format-plot t "set logscale xy")
(vgplot:plot (nreverse xs) (nreverse ys) "with linespoints")))
(plot-pi 1000000)Code language: Lisp (lisp)これだけで、πの有理近似の全体像が見えます。
- 355/113は5世紀の中国の数学者・天文学者である祖冲之（そ・ちゅうし、429–500）が発見したとされています。中国語では「密率（みつりつ）」と呼ばれ、4桁以下の分母を持つ有理数の中でπへの最良近似です。ヨーロッパで同じ値が独立に発見されたのは1585年のことで、祖冲之の発見から約1000年後でした。 – Milü – Wikipedia
- πが有理数ではない、つまり分数で表せないことを最初に厳密に証明したのはJohann Heinrich Lambert（1728–1777）で、1768年のことです。Lambertはtan(x)の連分数展開を使い、「xが0でない有理数ならtan(x)は有理数にならない」と示しました。tan(π/4)=1であることから、π/4は有理数ではありえないと結論しています。 – Proof that π is irrational – Wikipedia
- 有理数が実数直線上に稠密（dense）であるとは、任意の実数の近くに有理数が必ず存在することを意味します。この性質はディリクレの近似定理（後述）の前提にもなっており、「最も近い有理数」が存在しない理由でもあります。 – Rational approximation – AoPS Wiki
- Common Lispの
round関数は、最近接偶数への丸め（banker’s rounding）を行います。x*qが2つの整数のちょうど中間の場合、偶数の方に丸めます。πの近似では実際にそのような境界ケースはほぼ発生しませんが、一般的なアルゴリズムとして使うときは念頭に置いておく必要があります。 – Common Lisp HyperSpec: round - vgplotはVolker Sarodnick氏が開発したCommon Lispライブラリで、MATLABやOctaveのplot関数に似たAPIを提供します。quicklispから
(ql:quickload :vgplot)でインストールできます。別途gnuplotのインストールが必要です。 – volkers/vgplot – GitHub - 355/113は、分母が4桁以下の有理数の中でπへの最良近似です。誤差は約2.67×10⁻⁷で、πとの差はπの値の約0.0000085%にすぎません。次に分母が小さい競合候補（103993/33102）は分母が5桁になります。 – Milü – Wikipedia
- 連分数の理論は17〜18世紀にオイラー（Euler）やラグランジュ（Lagrange）らによって整備されました。収束分数が最良近似を与えるという定理（最良近似定理）は、ルジャンドル（Legendre）が1798年に証明しています。ある分数 p/q が収束分数であることの必要十分条件は |α – p/q| < 1/(2q²) を満たすことです。 – Simple continued fraction – Wikipedia
- 連分数の項に大きな数が現れると、その直前の収束分数が特によい近似になります。これは、次の項が大きいほど「次の補正量 1/(292+…)」が小さくなるためです。292はπの連分数に現れる係数の中で記録的に大きく、Wolfram MathWorldでもこの点が特記されています。 – Pi Continued Fraction – Wolfram MathWorld
- これはディリクレの近似定理（Dirichlet’s approximation theorem）から導かれる結論です。1842年にPeter Gustav Lejeune Dirichletが鳩の巣原理を使って証明しました。任意の無理数αに対して、|α – p/q| < 1/q² を満たす有理数p/qが無限に存在することを保証します。グラフの傾きが-2に近くなるのはこの定理の直接的な反映です。 – Dirichlet’s approximation theorem – Wikipedia