はじめに
アルキメデスは円に外接する正６角形、正12角形、正24角形、正48角形、正96角形と順に周囲の長さを計算することで、円周率の近似値を求めたとのことです。
tsujimotterさんノートブック "アルキメデスと円周率"
当時、小数点は未だ発明されておらず、整数の比で表現していたとのことです。
こんな制限があってもアルキメデスは、3と1/7より小さく、3と10/71より大きいことを求めていました。（小数点で表現すると、3.1428571428571429と3.1408450704225352の間）
アルキメデスの手法をそのまま進めると何所まで正確に円周率を計算できるかプログラミングの力を借りてをやってみました。
アルキメデスの手法
アルキメデスの手法そのものは、先のtujimotterさんのノートブックに正確に記載されていますので、興味のある方はそちらをご覧ください。ここでは、プログラム作りにあたって、必要な証明のみ記載します。
正n角形の周囲の長さをL、外接する円の半径を1とすると、円周は２πなので、正n角形を元にした円周率の近似値は L÷２となります。
【正6角形の場合】
L＝F(0)B×2×6
F(0)B ＝1/√3 ※ ∠F(0)OBが30°なので（△OF(0)Cは正三角形）
正6角形を元にした円周率の近似値は
L÷2＝1/√3×2×6÷2＝2×√3 ≒ 3.464 となる。
【正12角形の場合】
L＝F(1)B×2×12
F(1)B ＝ OB× F(0)B ÷（OB＋F(0)O）= F(0)B ÷（1＋√(1＋F(0)B× F(0)B))
正12角形をもとにした円周率の近似値は、
L÷2＝F(1)B×2×12÷2 ≒ 3.215 となる。
証明
F(1)Oと平行なF(0)を通る直線を引き、OBの延長線との交点をＡとする。
△OBF(1) と △ABF(0) は相似であるから、以下が成り立つ。
F(1)B：OB = F(0)B：（OB＋OA）
更に、△OF(0)Aは二等辺三角形であることから、F(0)O ＝ OA
F(0)Oはピタゴラスの定理を使って、√(OB×OB＋ F(0)B× F(0)B) より求められる。
【正ｎ角形の場合】
同様にして、正n角形を元にした円周率の近似値は以下のようになる。
L÷2＝F(n)B×2×n÷2 ＝ F(n-1)B÷(1＋√(1＋F(n-1)B× F(n-1)B))×ｎ
Pythonによる実装
以上の結果をPythonで実装すると下記のようになる。
import math
import itertools
def f(n):
if n == 0: return math.sqrt(3)/3
return f(n-1)/(1+math.sqrt(1+f(n-1)**2))
N = int(input())
for n in range(N):
print(f(n)*6*2**n)
ただし、正ｎ角形を正確に入力するのはメンドウなので、入力値Nは以下のようにしている。
0 … 正６角形を計算
1 … 正12角形まで計算
2 … 正24角形まで計算
以下同様。
上記のアルゴリズムではＯ(2**n）の計算量が必要になるため、次のようにDPで書き直してみました。また、２進数では小数点以下を正確に表現できない（※）ので、15桁ほど大きな数字にして、出力の際に整形するプログラムとしました。
※プログラミングで小数はヤバイんだぜ？という話。
import math
import itertools
def pi2s(pi):
return str(pi)[0] + "." + str(pi)[1:]
def regN(n):
return format(n, '10d') + ": "
R = 1732050807568877
B = 1000000000000000
pi = "3.14159265358979323846264338327950288"
N = int(input())
F = [0.0 for n in range(N)]
F[0] = B
for n in range(1, N):
F[n] = F[n-1]*R/(R+math.sqrt(F[n-1]**2+R**2))
for n, v in enumerate(F):
print(regN(6*2**n) + pi2s(int(v*6*2**n*B/R)))
print(format("pi: ", '>12') + pi[:17])
正3145728角形まで計算した結果です。
最後の行にπを出力しています。結構正確な値が計算できていることが分かります。
さいごに
アルキメデスの手法をプログラミングでなぞってみました。勝手な予想に反して結構早い段階(48角形)で既に3.14に近似できていて驚きでした。
これを、手でしかも整数の比（分数）で計算したアルキメデスには、ひたすら脱帽です。