祖沖之と円周率

この記事は以下の動画の補足記事です。

5 世紀ごろの中国で活躍した天文学者兼数学者、祖沖之は円周率について以下の結果を残したとされています。

  1. 約率と呼ばれる近似 $$\pi \fallingdotseq {\dfrac{22}{7}} \ \ (\ = 3.14285…)$$
  2. 密率と呼ばれる近似 $$\pi \fallingdotseq {\dfrac{355}{113}} \ \ (\ = 3.141592920…)$$
  3. 不等式 $$3.1415926 < \pi < 3.1415927$$

円周率 $\pi = 3.1415926535….$ なので小数点以下 $7$ 桁の精度であり、これは当時の世界での最高記録だったとされています。またその後約 900 年間、14 世紀ごろにインドの数学者マーダヴァが小数点以下 11 桁目まで求めるまでこの記録は破られなかったとされています。

しかし残念ながら、どのような計算を行ったのか、文献が残っていません。中国の歴史書、『隋書』律暦志よると、祖沖之は綴術という数学書を残していたようですが、彼の後を継いだもの達が理解できず、難しすぎて誰も学ばなくなり、次第に失われてしまったようです。冒頭の円周率の結果は『隋書』律暦志に残されていたものです。([上野])

残された情報はわずかではありますが、当時の暦の計算方法と合わせて、祖沖之の円周率の計算方法が推測されています ([曲], [杉本])。

暦と分数近似

現代の暦

暦とは、時間の流れを年・月・週・日といった単位に当てはめて数えるように体系付けたものです ([wiki1])。現在広く使われている暦は太陽暦といって、季節の周期 (春分・夏至・秋分・冬至の周期) を基本として 1 年を 365 日と定め (太陽年 [wiki2])、それを 12 の月に分けています (公転周期とは微妙に異なる)。

この周期は約 365.2422 日なので、これだと

$$(365.2422 -365) \times 4 = 0.9688$$

と、4 年で約 1 日足りなくなってしまいます。このずれを、4 年に 1 度 2月29日に閏日を設けることで調整しています。閏日がある年を閏年といいます。

しかしそれだと

$$\left(365.2422 -\left(365 + \frac{1}{4}\right) \right) \times 4= -0.0312$$

なので、4 年に約 0.03 日、400 年に約 3 日多く進んでしまいます。そこで 400 年に 3 回、100 の倍数かつ 400 の倍数でない年の閏年をスキップします。そうすると、400 年の間に $100 -3$ 回閏年が来るので

$$\begin{align*} & \left(365.2422 -\left(365 + \frac{1}{4} -\frac{3}{400}\right) \right) \times 400 \\ \ = & \left(365.2422 -\left(365 + \frac{97}{400}\right) \right) \times 400 \\ \ = & -0.12 \end{align*}$$

と、400 年に約 0.12 日多く進むことになります。これであれば 1 日ずれるのに 3000 年以上かかるので、これで良いだろうということになっているようです。

(ここで説明したのはグレゴリオ暦と呼ばれています。より誤差の少ない修正ユリウス暦というものもあります。)

ちなみにこの定め方は近似の効率という観点からではあまり効率が良くありません。例えば

\begin{align*} 365 + \frac{8}{33} &= 365.2424… \\ 365.2422 -365 + \frac{8}{33} &\fallingdotseq -0.00022 \end{align*}

なので、33 年に 8 回閏年を設ければ、約 4500 年に 1 日しかずれなくなります。他にも

\begin{align*} 365 + \frac{31}{128} &= 365.2421875 \\ 365.2422 -365 + \frac{31}{128} &= 0.0000125 = \frac{1}{80000} \end{align*}

なので、128 日に 31 回閏年を設ければ、例えば基本的に 4 年に 1 回閏年を設け 128 年に一度は設けないようにすれば、8 万年に 1 日しかずれなくなります。

古代中国の暦

古代中国ではこれとは異なり、太陰太陽暦と呼ばれるタイプの暦を用いていました。まず月の満ち欠けの周期(約 29.5306 日)を基準としてひと月を 29 日または 30 日と定め、29 日の月と 30 日の月を交互に 6 回繰り返したものを 1 年と定めます。

しかしこれだと 1 年の日数は $6 \times 59 = 354$ なので、1 年間で 365 日から約 11 日ずれてしまいます。太陰太陽暦ではこのずれを、数年に一度閏月を設けることで修正します。

ではどのくらいの頻度で閏月を設ければ良いのでしょうか。上記のように計算しても良いですが、古代中国ではもっと効率の良い近似ができる方法で計算していました。

閏月の計算と分数近似

季節の周期を $T$ $(\fallingdotseq 365.2422)$ 日、月の満ち欠けの周期を $B$ $(\fallingdotseq 29.5306)$ 日とします。 (これらの値は現代の計算結果であって、当時の中国で求められていた値とは異なる可能性があります。[曲])

閏月をどのくらいの頻度で設ければよいのかという問題は、$\dfrac{T}{B}$ を分数で近似するという問題と考えられます。

$$\frac{T}{B} = \frac{365.2422}{29.5306} \fallingdotseq 12.3683$$

なのでおおよそ $12$ ですが、残りの $0.368…$ を $\dfrac{1}{3}$ と近似することは、3 年に一度閏月を設けることを意味します。そのとき、季節の周期のずれは

$$\frac{T}{B} -\left(12 +\dfrac{1}{3} \right)$$

となります。残りの部分を

$$\dfrac{7}{19} \fallingdotseq 0.368421$$

と近似することは、19 年に 7 度閏月を設けることを意味します。このとき、季節の周期のずれは

$$\frac{T}{B} -\left(12 +\dfrac{7}{19} \right)$$

となります。

ではこのような分数をどのように求めれば良いでしょうか?グレゴリオの暦のように求めても良いですが、それだと $\dfrac{8}{33}$ や $\dfrac{31}{128}$ のような、より近似精度が高く、より “単純に見える” 分数を見落としてしまいます。

古代中国ではそのような分数を見つける方法を知っていました。

調日法

分数の見つけ方の説明の前に、”単純に見える” 分数について補足しておきます。”単純に見える” とは単に分母が小さいという意味です。古代中国では、実数 $\alpha$ に対して、最も $\alpha$ を精度良く近似する分数を分母が小さいものから順に見つける方法を知っていました。

それは何承天 (370年 – 447年) が編み出した、調日法と呼ばれる手法です。調日法とは、分数が

$$\frac{a}{b} < \frac{c}{d}$$

を満たすときに

$$\frac{a}{b} < \frac{a + c}{b + d} < \frac{c}{d}$$

が成り立つ、という性質を利用して、実数 $\alpha$ に対して

$$\frac{a}{b} < \alpha < \frac{c}{d}$$

という不等式が得られた時に $\alpha$ と $\dfrac{a + c}{b + d}$ の大小を比較して、より精度の高い分数近似を得る方法です。

例えば $\dfrac{T}{B} -12$ の近似分数を求めるには以下のようにします。まず不等式

$$\frac{0}{1} < \frac{T}{B} -12 < \frac{1}{1}$$

から始めて、$\dfrac{0+1}{1+1} = \dfrac{1}{2}$ と $\dfrac{T}{B} -12$ を比較すると不等式

$$\frac{0}{1} < \frac{T}{B} -12 < \frac{1}{2}$$

が得られます。もう一度繰り返して、$\dfrac{0 + 1}{1 +2} = \dfrac{1}{3}$ と $\dfrac{T}{B} -12$ を比較すると不等式

$$\frac{1}{3} < \frac{T}{B} -12 < \frac{1}{2}$$

が得られます。これを繰り返すと区間がどんどん短くなっていくので、$\dfrac{T}{B} -12$ に近い方の端点を取れば $\dfrac{T}{B} -12$ の近似分数が得られます。

このようにして近似分数を求められますが、

  1. 本当に不等式 $$\frac{a}{b} < \frac{a + c}{b + d} < \frac{c}{d}$$ が成り立つのか
  2. 本当に “分母の小さい” 近似分数を見つけられるのか

をこれ以降で確認していきます。

(補足) 調日法は単に近似精度を上げる目的のみにも有効なので、当時の中国において “分母の小さい” 分数近似が求められることを認識していたのか、定かではありません。初等的に示すことができるので、当時の中国の数学の水準的に証明することは可能だったと思いますが、詳しく調査できていません。

実数の分数近似

前置き

分子分母が整数である分数 $\dfrac{a}{b}$ が既約であるとは、$a$ と $b$ が共通因子を持たないこと、つまり $a$ と $b$ が互いに素であることとします。ただし $a = 0$ の場合は $b = 1$ のときに既約であるとします。

以下この記事では、有理数を $\dfrac{a}{b}$ と表したとき、特に明記しない限り既約であり、かつ $b > 0$ であるとします。

分数表示の一意性

有理数を分数で表示するとき、その方法は一通りではありません。例えば

$$\frac{1}{2} = \frac{2}{4} = \frac{-2}{-4}$$

のように、値としては等くても表し方は複数あります。ただし、既約かつ分母が正である、という条件のもとで表示は一意的定まることが以下のようにわかります。

既約分数 $\dfrac{p_1}{q_1}$, $\dfrac{p_2}{q_2}$ が有理数として等しいとします。このとき

$$0 = \dfrac{p_1}{q_1} -\dfrac{p_2}{q_2} = \frac{p_1q_2 -p_2q_1}{q_1q_2}$$

なので、$p_1q_2 = p_2q_1$ となります。よって $p_1$ は $p_2q_1$ を割り切ります ($\dfrac{p_2q_1}{p_1}$ が整数となる) が、$p_1$ と $q_1$ は互いに素なので $p_1$ は $p_2$ を割り切ります。同様に $p_2$ が $p_1q_2$ を割り切ること、$p_2$ と $q_2$ が互いに素であることから $p_2$ は $p_1$ を割り切ります。したがって $p_1 = p_2$ または $p_1 = -p_2$ となります。

ここで、$p_1q_2 = p_2q_1$ と $q_1 > 0$, $q_2 > 0$ から $p_1$ と $p_2$ の符号は等しいです。したがって $p_1 = p_2$ となります。さらに $p_1q_2 = p_2q_1$ から $q_1 = q_2$ も従います。

良い近似とは

実数 $\alpha$ を分数 $\dfrac{a}{b}$ で近似するとき、どのような尺度でその良さを決めれば良いでしょうか。近似の精度の良さはもちろん重要ですが、分母を大きくすればいくらでも近似の精度をあげることができます。その分、分数としては複雑になってしまいます。

よって近似の精度がよく、より単純な分数が良い近似分数であると考えたいところですが、分数の単純さをどう定義すれば良いでしょうか?色々な考え方があると思いますが、ここでは分母が小さい分数の方が単純である、と考えることにします。

この考え方のもと、実数 $\alpha$ の近似分数 $\dfrac{a}{b}$ が良い近似分数であることを、$b \geq d$ を満たし $\dfrac{a}{b}$ と異なる任意の分数 $\dfrac{c}{d}$ に対して

$$\left|\alpha -\dfrac{a}{b}\right| < \left|\alpha -\dfrac{c}{d}\right|$$

を満たすものとします。つまり、$\dfrac{a}{b}$ より分母の小さい分数で、$\dfrac{a}{b}$ より近似精度が良いものが存在しないとき、$\dfrac{a}{b}$ を $\alpha$ の良い近似分数であるとします。(これは一般に第一種最良近似分数と呼ばれます。)

例えば $\dfrac{97}{400}$ は、より分母の小さい分数 $\dfrac{8}{33}$ が $0.2422$ をより精度よく近似するため、$0.2422$ の良い近似分数ではありません。自明ではありませんが、$\dfrac{8}{33}$ は良い近似分数です。この後で説明しますが、それは調日法によって得られることから従います。

中間分数と不等式

分数 $\dfrac{a}{b} < \dfrac{c}{d}$ に対して、不等式

$$\dfrac{a}{b} < \dfrac{a +c}{b + d} < \dfrac{c}{d}$$

が成り立ちます。右側の不等式は

$$\begin{align*} \frac{c}{d} -\frac{a +c}{b + d}&= \frac{c(b+d) -d(a+c)}{d(b+d)} \\ &= \frac{cb -ad}{d(b+d)} > 0 \end{align*}$$

から成り立ちます。ここで最後の不等式は

$$\dfrac{a}{b} < \dfrac{c}{d} \Rightarrow ad < cb$$

を用いました。左側も同様に計算すれば良いです。$\dfrac{a +c}{b + d}$ を中間分数と呼びます。

重要な補題とその証明

重要な補題をいくつか証明します。

補題 1.

2 つの分数 $\dfrac{a}{b}, \dfrac{c}{d}$ は $\dfrac{a}{b} < \dfrac{c}{d}$, $bc -ad = 1$ を満たすとする. このとき中間分数 $\dfrac{a +c}{b + d}$ は $\dfrac{a}{b} < \dfrac{x}{y} < \dfrac{c}{d}$ を満たす分数のうち分母が最小のものであり, 分母が最小のものは $\dfrac{a +c}{b + d}$ のみである.

(証明)

$b + d$ が分母の最小値であることを示します。$\dfrac{a}{b} < \dfrac{x}{y}$ の両辺に $yb$ をかけて整理すると $xb -ay > 0$ となります。左辺は整数なので、$xb -ay \geq 1$ となります。$\dfrac{x}{y} < \dfrac{c}{d}$ も同様にして $cy -xd \geq 1$ となります。ここで

$$\begin{align*} v &= xb -ay \\ u &= cy -xd \end{align*}$$

とおくと、$cb -ad = 1$ が成り立つことから $x, y$ について解いて

$$\begin{align*} x &= au + cv \\ y &= bu + dv \end{align*} $$

となります。ここで、$u \geq 1, v \geq 1$ から $y \geq b +d$ でなければなりません。よって $b + d$ が分母の最小値となります。

最後に分母が最小のものが $\dfrac{a +c}{b + d}$ ただ一つであることを示します。分母が最小という条件から $y = b + d$ となるので

$$\begin{align*} & b + d = bu + dv \\ \Leftrightarrow \ \ & b(u -1) + d(v -1) = 0 \end{align*}$$

となりますが、$b, d, u, v \geq 1$ なので $u = v = 1$ となります。よって

$$x = au + cv = a + c$$

となり、分母が最小のものが $\dfrac{a +c}{b + d}$ ただ一つであることがわかります。$\Box$

分母が最小であることから、$\dfrac{a +c}{b + d}$ が既約であることもわかります。(そうでなければ約分をしたもののほうが分母が小さい。)

ちなみに $ad -bc = 1$ から $\dfrac{a}{b}$, $\dfrac{c}{d}$ は既約であることが従います。実際、もし $\dfrac{a}{b}$ が既約でないなら、$a$, $b$ はある共通因子 $p$ をもち、$ad -bc$ は $p$ の倍数になってしまいます。$\dfrac{c}{d}$ についても同様です。

補足

線形代数の知識があると証明の見通しが良くなります。まず、$\mathbb{R}^2$ の格子点 $(x, y)$ のうち、$y \geq 1$ かつ $x$ と $y$ が互いに素であるものと、有理数が 1 対 1 に対応することに注意します。(互いに素であることと、原点と $(x, y)$ を結んだ線分の間に格子点が存在しないことが同値です。)

そこで、行列

$$\begin{align*} \begin{pmatrix}c & a \\ d & b\end{pmatrix} \end{align*}$$

を考えると、これは (左からかけることで) ${}^t (1, 0)$ を ${}^t (c, d)$ に移し、${}^t (0, 1)$ を ${}^t (a, b)$ に移します。さらに全ての成分が整数なので格子点を格子点に移し、逆行列

$$\begin{align*} \frac{1}{cb -ad}\begin{pmatrix}b & -a \\ -d & c\end{pmatrix} = \begin{pmatrix}b & -a \\ -d & c\end{pmatrix} \end{align*}$$

も全ての成分が整数なので、この行列によって格子点が 1 対 1 に対応します。特に、第一象限はベクトル ${}^t (c, d)$ と ${}^t (a, b)$ の間の領域 (ベクトルを伸ばしてできる領域) に移るので、$x \geq 1, y \geq 1$ を満たす格子点は ${}^t (c, d)$ と ${}^t (a, b)$ の間にある格子点に移ります。そしてこれは $\dfrac{a}{b} < \dfrac{x}{y} < \dfrac{c}{d}$ を満たす有理数 $\dfrac{x}{y}$ と 1 対 1 に対応します。

この対応によって格子点 $(u, v)$ が有理数 $\dfrac{cu + av}{du + bv}$ に対応し、$u \geq 1, v \geq 1$ を満たす格子点のうち、$u$ と $v$ が互いに素であるものと、$\dfrac{a}{b}$ と $\dfrac{c}{d}$ の間の分数が 1 対 1 に対応します。これから補題の主張が直ちに従います。ついでに、$\dfrac{a}{b}$ と $\dfrac{c}{d}$ が既約なら $\dfrac{c + a}{d + b}$ も既約であることがわかります。

ファレイ数列と中間分数

自然数 $n$ に対して、$0$ 以上 $1$ 以下の分数で、分母が $n$ 以下のものを小さいものから順に並べたものを $F_n$ と表します。$F_n$ をファレイ数列と言います。例えば

$$F_4 = \left\{ \frac{0}{1}, \frac{1}{4}, \frac{1}{3}, \frac{1}{2}, \frac{2}{3}, \frac{3}{4}, \frac{1}{1} \right\}$$

です。$0$ 以上 $1$ 以下という条件を外し、分母が $n$ 以下の全ての既約分数を小さいものから順に並べたものを $\hat{F}_n$ と表します。

ファレイ数列を用いると $bc -ad = 1$ という条件を図形的に意味づけることができます。

補題 2.

既約分数 $\dfrac{a}{b} < \dfrac{c}{d}$ に対して, $bc -ad = 1$ を満たすことと, $\hat{F}_{\max\{b, d\}}$ 内で隣合うことは同値である.

(証明) ($\Rightarrow$) $bc -ad = 1$ を満たすとします。補題 1 から、$\dfrac{a}{b} < \dfrac{x}{y} < \dfrac{c}{d}$ を満たす分数 $\dfrac{x}{y}$ のうち分母が最小のものは $\dfrac{a+c}{b+d}$ で、かつ $\max\{b, d\} < b + d$ なので $\dfrac{a}{b}$ と $\dfrac{c}{d}$ の間には $\hat{F}_{\max\{b, d\}}$ の要素が存在しません。

($\Leftarrow$) ${\max\{b, d\}}$ に関する帰納法で証明します。${\max\{b, d\}} = 1$ のとき、$\hat{F}_1$ の要素は $\dfrac{n}{1}$ という分数のみで、$(n+1) \cdot 1 -n \cdot 1 = 1$ なので成立します。${\max\{b, d\}} = n$ まで成立すると仮定し、${\max\{b, d\}} = n + 1$ のときにも成立することを示します。

まず、$\hat{F}_{n+1}$ に含まれて $\hat{F}_n$ に含まれない分数の集合 $\hat{F}_{n+1} \setminus \hat{F}_n$ の要素同士は $\hat{F}_{n+1}$ で隣り合いません。それは以下のようにしてわかります。

$\hat{F}_{n+1} \setminus \hat{F}_n$ の要素は $\dfrac{k}{n+1}$ という形の既約分数です。これが $\hat{F}_{n+1} \setminus \hat{F}_n$ の要素と隣合うとすれば、それは $\dfrac{k-1}{n+1}$ または $\dfrac{k+1}{n+1}$ です。どちらの場合も同じなので、

$$\dfrac{k}{n+1}, \dfrac{k+1}{n+1} \in \hat{F}_{n+1} \setminus \hat{F}_n$$

のときに、$\dfrac{k}{n+1}$ が $\dfrac{k+1}{n+1}$ と隣合うのかを考えます。

$\dfrac{k}{n+1} \in \hat{F}_{n+1} \setminus \hat{F}_n$ から、整数 $m$ と $1 \leq k^{\prime} \leq n$ が存在して $\dfrac{k}{n+1} = m + \dfrac{k^{\prime}}{n+1}$ と表すことができます。このとき

$$\dfrac{k}{n+1} = m +\dfrac{k^{\prime}}{n+1} < m +\frac{k^{\prime}}{n} < m +\frac{k^{\prime}+1}{n+1} = \dfrac{k+1}{n+1} $$

のように間に $m +\dfrac{k}{n} \in \hat{F}_n$ が挟まるため、$\dfrac{k}{n+1}$ と $\dfrac{k+1}{n+1}$ は隣り合いません。実際

\begin{align} \frac{k^{\prime}}{n} -\dfrac{k^{\prime}}{n+1} &= \frac{(n+1)k^{\prime} -nk^{\prime}}{n(n+1)} = \frac{k^{\prime}}{n(n+1)} > 0 \\ \frac{k^{\prime}+1}{n+1} -\frac{k^{\prime}}{n} &= \frac{n(k^{\prime}+1)-(n+1)k^{\prime}}{n(n+1)} = \frac{n-k^{\prime}}{n(n+1)} \geq 0 \\ \end{align}

であり、$n-k^{\prime} = 0$ ならば $\dfrac{k+1}{n+1}$ は整数なので仮定に反します。

よって $\dfrac{k}{n+1} \in \hat{F}_{n+1} \setminus \hat{F}_n$ の要素の両隣の分数は $\hat{F}_n$ の要素です。それらの分数を $\dfrac{a}{b} < \dfrac{c}{d}$ とすると、これらは $F_n$ で隣合うので、仮定から $bc -ad = 1$ で、補題 1 からその中間分数 $\dfrac{a+c}{b+d}$ は $\dfrac{k}{n+1}$ に一致します。

ここで、$\hat{F}_{n+1}$ で隣合う分数 $\dfrac{a}{b}$, $\dfrac{a+c}{b+d}$, $\dfrac{c}{d}$ は

\begin{gather*}b(a+c) -a(b+d) = bc -ad = 1\\ c(b+d) -d(a+c) = bc -ad = 1\end{gather*}

を満たすので、$n+1$ でも条件を満たすことが示されました。$\Box$

小数部分のみの調日法

調日法は、最初の不等式をどう与えるかによりますが、最初の不等式を整数 $n$ による不等式 $n < \alpha < n+1$ から始める限り、$\alpha$ の小数部分 $\alpha -\lfloor \alpha \rfloor$ のみを対象にしても同じ結果が得られます。実際、$n + \dfrac{a}{b} = \dfrac{nb + a}{b}$ と $n + \dfrac{c}{d} = \dfrac{nd + c}{c}$ の中間分数は

$$\dfrac{nb + a + nd + c}{b + c} = n + \frac{a + c}{b + d}$$

となります。つまり、整数 $n$ 平行移動してから中間分数をとっても、中間分数をとってから整数 $n$ 平行移動しても同じ結果が得られます。

よって $0 \leq \alpha < 1$ を満たす実数に対して議論を展開すれば十分です。

良い近似分数は調日法で計算できる

調日法のように中間分数をとる操作を繰り返すことで、良い近似分数の列を求められます。

定義から、実数 $\alpha$ が不等式 $\dfrac{a}{b} < \alpha < \dfrac{c}{d}$ を満たし、両端の分数が $\hat{F}_{\max\{b, d\}}$ で隣合うならば、両端の分数のうち、$\alpha$ に近いほうが良い分数近似です。

補題 1、補題 2 から、中間分数 $\dfrac{a + c}{b + d}$ と $\dfrac{a}{b}$, $\dfrac{c}{d}$ は $\hat{F}_{b + d}$ で隣合います。区間 $\left(\dfrac{a}{b}, \dfrac{c}{d}\right)$ を中間分数 $\dfrac{a + c}{b + d}$ で 2 つに分け、$\alpha$ が含まれる方の区間を $\left(\dfrac{a^{\prime}}{b^{\prime}}, \dfrac{c^{\prime}}{d^{\prime}}\right)$ とおけば、両端の分数は $\hat{F}_{b + d}$ で隣合うので、$\alpha$ に近いほうが良い近似分数になります。($\alpha = \dfrac{a + c}{b + d}$ となる場合もありますが、その場合はその時点で調日法が終了します。)

調日法の最初の不等式を、$\alpha$ に対して整数 $n$ で、$n < \alpha < n+1$ を満たすもので取ります。これは $\hat{F}_1$ で隣合い、$n$ と $n + 1$ のうち $\alpha$ に近いほうが良い分数近似です。よって上記の方法を繰り返して得られる分数は全て良い分数近似です。

円周率の分数近似

劉徽は以下の結果を残したとされています。

\begin{gather*}\pi \fallingdotseq 3.1416 \\ 3.141024 < \pi < 3.142704\end{gather*}

1 つ目を劉徽の近似値、2 つ目を劉徽の不等式と呼ぶことにします。祖沖之がこれを知っていたことを前提に、円周率を分数で近似しましょう。

円周率 $\pi$ の値は正確にはわからないので、代わりに近似値 $\alpha = 3.1416$ を分数近似します。

$0 < \alpha -3 < 1$ から始めて、中間分数

$$\dfrac{0 + 1}{1 + 1} = \dfrac{1}{2} = 0.5$$

が得られますが、劉徽の結果から $0 < \alpha -3 < \dfrac{1}{2}$ がわかります。次の中間分数は

$$\frac{0 + 1}{1 + 2} = \frac{1}{3} = 0.333…$$

となります。再び劉徽の結果から $0 < \alpha -3 < \dfrac{1}{3}$ がわかります。これを繰り返して中間分数を求めていっても良いですが、次のようにすると効率が良いです。

命題.

$\dfrac{a}{b} < \dfrac{c}{d}$ とする. 任意の実数 $0 \leq y < x$ に対して

$$\frac{a + c y}{b + d y} < \frac{a + c x}{b + d x}$$

が成り立つ.

(証明) $b > 0$, $d > 0$ に注意して計算すれば良い。

\begin{align*} &\frac{a + c x}{b + d x} -\frac{a + c y}{b + d y} \\ = \ & \frac{(a + c x)(b+dy) -(a + c y)(b + d x)}{(b + d x)(b + dy)} \\ = \ & \frac{cbx + ady -cby -adx}{(b + d x)(b + dy)} \\ = \ & \frac{(cb -ad)(x-y)}{(b + d x)(b + dy)} > 0\\ \end{align*}

$\Box$

$\dfrac{a + c x}{b + d x}$ は $x \to \infty$ で $\dfrac{c}{d}$ に収束するので、$x > 0$ で

$$\frac{a}{b} < \frac{a + c x}{b + d x} < \frac{c}{d}$$

であり、$x$ が小さいほど $\dfrac{a}{b}$ に近く、$x$ が大きいほど $\dfrac{c}{d}$ に近くなります。

$x$ を変数として、$\dfrac{0}{1}$ と $\dfrac{1}{1}$ の中間分数 $\dfrac{0 x + 1}{x + 1}$ を考えます。この分数は、$x$ が小さいほど $1$ に近く、$x$ が大きいほど $0$ に近づきます。そして $\dfrac{1}{x + 1} < \alpha -3$ を満たす最小の正の整数 $x$ を求めれば

$$\dfrac{1}{x + 1} < \alpha -3 < \dfrac{0 +1}{(x-1) + 1} = \dfrac{1}{x}$$

という不等式が得られます。

約率を求める

まず、

$$\dfrac{1}{x + 1} = \alpha -3 = 0.1416$$

という方程式を解くと $x \fallingdotseq 6.06214$ が得られます。よって、$x = 7$ あたりが条件を満たすものだと予想できます。実際に計算すると

$$0.125 = \frac{1}{8} < \alpha -3 < \frac{1}{7} \fallingdotseq 0.142857$$

となります。劉徽の不等式から

$$\frac{1}{8} < \pi -3 < \frac{1}{7}$$

も成り立ち、$\dfrac{1}{8}$ のよりも $\dfrac{1}{7}$ の方が $\pi -3$ に近いこともわかるので、$\dfrac{1}{7}$ は $\pi-3$ の良い近似分数であることがわかります。

密率を求める

次は

$$\frac{1}{8} < \alpha -3 < \frac{1}{7}$$

に対して同様に、$\dfrac{1 + x}{8 + 7x}$ という分数を考えましょう。これは $x$ が小さいほど $\dfrac{1}{8}$ に近く、$x$ が大きいほど $\dfrac{1}{7}$ に近くなります。

$\dfrac{1 + x}{8 + 7x} = 0.1416$ を解いてみましょう。すると $x \fallingdotseq 15.09$ となります。$x = 15$ だとこの分数は $0.1416$ より小さく、$x = 16$ だと $0.1416$ より大きくなります。実際に計算すると

$$0.1415929 \fallingdotseq \frac{16}{113} < \frac{17}{120} \fallingdotseq 0.1416667$$

となります。ここで $3 + \dfrac{16}{113} = \dfrac{355}{113}$ なので、祖沖之が求めた密率が現れます。

しかし上記で求めたものは $0.1416$ の近似分数なので、実際の円周率とどのくらい近いのかは分かりません。もう一つ前の近似分数は

$$\frac{15}{106} = 0.141509$$

であることを考えると、$3 + \dfrac{15}{106}$, $3 + \dfrac{16}{113}$, $3 + \dfrac{17}{120}$ のなかで $3 + \dfrac{16}{113}$ が最も円周率に近そうですが、本当にそうなのでしょうか。また、どのくらい精度良く近似できているのでしょうか。

祖沖之の求めた不等式

最後に、祖沖之が以下の不等式をなぜ、どのように求めたのかについて考察します。

$$3.1415926 < \pi < 3.1415927$$

求めた理由はおそらく、区間 $\left(3 + \dfrac{15}{106}, 3 + \dfrac{16}{113}\right)$ に含まれるのか、$\left(3 + \dfrac{16}{113}, 3 + \dfrac{17}{120}\right)$ に含まれるのかを確認したかったのだと思います。調日法を進めるにはそれを知る必要があります。劉徽の結果から、$3 + \dfrac{16}{113}$ が最も近いだろうことは予測がついていたと思いますが、どちらに含まれるかは微妙な問題です。

ではどのように不等式を求めたのでしょうか。情報がないため、劉徽の方法で求めたと仮定します。劉徽の方法は円を正多角形の面積で近似する手法ですが、平方根の計算が含まれるため、平方根の値をある桁数で打ち切る必要があります。平方根の 15 桁以降を切り捨て、正 6 角形から始めて正 24576 角形まで計算すると、以下の不等式が得られます。

$$3.1415926017… < π < 3.1415926918…$$

ちなみに平方根 15 桁以降を切り捨てた影響で、計算誤差が溜まっていき、正 49152 角形では以下のような不等式が得られます。

$$3.1415924815… < π < 3.1415923614…$$

特に、単調増加するはずの左辺が減少してしまいます。

ちなみにもっと平方根の計算精度を上げると、以下のような不等式が得られます。

\begin{align*} 3.1415926188… < &π < 3.1415927212… \quad (24576 \textrm{ 角形}) \\ 3.1415926394… < &π < 3.1415926600… \quad (49152 \textrm{ 角形})\end{align*}

よって、劉徽の方法で計算したなら、平方根を 14 桁以上の精度で計算し、少なくとも正 24576 角形か正 49152 角形まで面積を計算したのではないかと考えられます。

さらなる分数近似

再掲しますが、

$$3 + \dfrac{16}{113} = 3.141592920…$$

です。よって祖沖之の不等式から

$$\pi \in \left( 3 + \dfrac{15}{106}, 3 + \dfrac{16}{113} \right)$$

であることがわかります。この中間分数は

$$3 + \frac{31}{219} = 3.1415525…$$

となります。これは $3 + \dfrac{16}{113}$ に比べて円周率からかなり遠いです。調日法を続けていくと、区間の左側がどんどん狭まっていきますが、それが祖沖之の不等式の左辺 $3.1415926$ を超える分母が最も小さい分数は

$$\frac{86953}{27678} = 1415926006214323…$$

です。祖沖之の不等式から、これが $3 + \dfrac{16}{113} = \dfrac{355}{113}$ よりも精度の良い近似であることがわかりますが、かなり複雑です。

計算結果を github で公開しています。気になる方はリンクから参照ください。

参考文献

[曲] 曲 安京 (著), 城地 茂 (訳). 祖沖之は、 如何に円周率 π = 355/113 を得たか ?

[杉本] 杉本 敏夫. 中国古代の周率 (下)

[上野] 上野 健爾. 円周率が歩んだ道

[H] G. H. HARDY, E. M. WRIGHT. AN INTRODUCTION TO THE THEORY OF NUMBERS

[Wiki] Wikipedia.

[Wiki2] Wikipedia. 太陽年

[高木] 高木貞治 初等整数論講義第2版 (Line Segments) §22 近似分数の特徴

[橋本] 橋本竜太. 実数の有理数近似と連分数展開

シェアする