뉴턴의 이항정리를 활용한 원주율 계산법

1. 개요

뉴턴 이전의 사람들은 원주율을 계산하는 방식은 아르키메데스가 만든 방식을 사용했다. 그의 방식은 반지름이 1인 원에 내접하는 정N면체와 외접하는 정X면체의 넓이를 구해서 원주율값을 구해나가는 것이었다. 예를들어 원에 외접하는 정사면체의 넓이는 4이고 내접하는 정사면체의 넓이는 2이므로 원주율은 2와 4 사이에 존재하는 것은 쉽게 알 수 있다. 이를 통해 아르키메데스는 정 96면체를 작도해서 3.16정도의 근사값을 얻었다. 뤼돌프 판 쾰런은 정 2^26면체를 활용해서 원주율의 35자리까지 근사했다. 이 때 걸린 시간은 28년이라고 한다.

 

그러다가 뉴턴은 이항정리를 활용해 원주율 계산방식을 매우 효율적으로 정확하게 계산하게 되면서 아무도 사용하지 않게 되었다. 지금은 라이프니츠의 교대급수를 사용하지만 결과적으로 모두 연결된 공식이다.

 

2. 코드

파이썬으로 해당 부분을 작성했고 이를 몬테카를로 시뮬레이션으로 계산하는 원주율과 비교해봤다.

 

#몬테카를로 시뮬레이션
import random
import matplotlib.pyplot as plt
import time
import numpy as np
def calc_time_pi_by_simulation(n=100000):
    start=time.time() # starting time
    count=0
    for k in range(n):
        x=random.random() # 0~1사이의 값 반환
        y=random.random()
        if x**2+y**2<=1:
            count+=1
    end=time.time() # end time
    approximated_pi=count/n*4
    return approximated_pi, end-start


calc_time_pi_by_simulation(n=100000000)

빨간색이 원의 1/4에 들어오는 점을 표시하고 나머지는 바깥 쪽에 있는 것이다.

위의 코드를 시행하면 컴퓨터에 따라 다르겠지만 나의 경우 40초가 걸렸다.

 

# 뉴턴의 이항정리를 활용한 원주율 계산 함수
import math
def calc_newton_coef(r=1/2,k=5):
    denominator=[math.factorial(i) for i in range(k)]
    numerators=[1,r]
    for n in range(2,k):
        numerators.append((numerators[n-1])*(r-n+1))
    coefs=[numerators[i]/denominator[i] for i in range(k)]
    return coefs


def calculate_pi(k=5):
	'''
    k= k-1 차항까지 계산에 포함할 것인지를 지정한다.
    '''
    start=time.time()
    newton_coefs=np.array(calc_newton_coef(r=1/2,k=k))
    weights=[(2*i+1) for i in range(0,k)] # 0, 2, 4, 6, 8, ...
    weights=[((1/2)**w)/w for w in weights]
    values=weights*newton_coefs
    values[1:]=[w if w<0 else -w for w in values[1:]]
    pi=(sum(values)-(3**(1/2))/8)*12
    end=time.time()
    return pi,end-start
    
approximated_pi=calculate_pi(k=100) # 파이썬에서 계산이 가능한 한계치.
# math.factorial부분을 좀 더 다듬으면 좀 더 빠르고 효과적일거 같은데 생략한다.

 

3. 걸리는 시간 비교

import seaborn as sns
def visualize_graph_for_pi(n=5):
    pi=math.pi
    result=[]
    for i in range(n):
        apx_pi,seconds=calc_time_pi_by_simulation(n=100000)
        error=abs(apx_pi-pi)
        result.append(['simulation',error,seconds])
        apx_pi,seconds=calculate_pi(k=5)
        error=abs(apx_pi-pi)
        result.append(['newton',error,seconds])
    result=pd.DataFrame(result,columns=['method','error','time(sec)'])
    sns.boxplot(data=result,x='method',y='time(sec)');plt.show()
    sns.boxplot(data=result,x='method',y='error');plt.show()
    sns.scatterplot(data=result,x='time(sec)',y='error',hue='method')
    plt.show()
    return result.groupby('method').mean()


visualize_graph_for_pi(n=5)

시간은 약 100배 가량 차이가 난다.

오차 역시 약 100배 가량 차이가 난다.

5회 반복을 통해 얻은 평균 오차값과 소요 시간

 

4. 결론

뉴턴 방식은 매우 빠르게 원주율값을 시뮬레이션에 비해 계산하며 4차항까지 계산해도 그 오차값은 1.7E-5밖에 안되므로 실생활에서 문제가 안될 수준에 이른다.