소개

아르키메데스는 원에 외접하는 정 6각형, 정 12각형, 정 24각형, 정 48각형, 정 96각형과 순서대로 주위의 길이를 계산하는 것으로, 원주율의 근사치를 구했다고 합니다.
tsujimotter씨 노트북 "아르키메데스와 원주율"
당시 소수점은 아직 발명되지 않았고 정수의 비율로 표현했다는 것입니다.
이런 제한이 있더라도 아르키메데스는 3과 1/7보다 ​​작고 3과 10/71보다 큰 것을 요구했습니다. (소수점으로 표현하면 3.1428571428571429와 3.1408450704225352 사이)
아르키메데스의 수법을 그대로 진행하면 몇 곳까지 정확하게 원주율을 계산할 수 있는지 프로그래밍의 힘을 빌려서 해 보았습니다.

아르키메데스 기법

아르키메데스의 방법 자체는 앞의 tujimotter씨의 노트북에 정확하게 기재되어 있으므로, 관심있는 분은 그쪽을 봐 주세요. 여기에서는 프로그램 만들기에 있어서 필요한 증명만 기재합니다.
정n각형의 주위의 길이를 L, 외접하는 원의 반경을 1로 하면, 원주는 2π이므로, 정n각형을 바탕으로 한 원주율의 근사치는 L÷2가 됩니다.

【정 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의 연장선과의 교점을 A로 한다.
△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각형의 경우】

마찬가지로, 양의 n각형에 기초한 원주율의 근사값은 다음과 같다.
L÷2=F(n)B×2×n÷2 = F(n-1)B÷(1+√(1+F(n-1)B×F(n-1)B))×n

파이썬으로 구현

이상의 결과를 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각형을 정확하게 입력하는 것은 멘도우이므로, 입력값 N은 이하와 같이 하고 있다.
0 … 정 6각형을 계산
1 … 정 12각형까지 계산
2 … 정24각형까지 계산
이하 동일.
위의 알고리즘에서는 O(2**n)의 계산량이 필요하기 때문에 다음과 같이 DP로 다시 작성해 보았습니다. 또, 2진수에서는 소수점 이하를 정확하게 표현할 수 없다(※) 때문에, 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에 근사할 수 있어 놀라웠습니다.
이것을, 손으로 게다가 정수의 비(분수)로 계산한 아르키메데스에는, 오로지 탈모입니다.