BigFloat로 Pi를 구해보자-11. Pi 구하기: Machin 급수

Posted by 적분 ∫2tdt=t²+c
2012.11.22 01:38 프로그래밍/Multi precision



저번 연재에서 아크탄젠트 급수를 이용해서 Pi를 계산해보았는데, 머리를 조금 더 굴리면 훨씬 더 빠른 알고리즘을 만들어낼 수 있습니다.

가장 대표적인 공식은 다음과 같아요.

John Machin이 1706년에 발견한 급수인데, Machin의 이름을 따서 Pi를 계산해내는 일련의 아크탄젠트 공식들을 Machin-like 공식이라고 부릅니다.[각주:1]


그냥 지나가면 재미없으니, 저 공식이 어떻게 유도되었는지 거슬러 올라가봅시다. (복소수를 이용한 증명을 소개하니, 이해하려면 복소수에 대한 기초적인 지식이 필요해요..)


...(식1)

복소수 지수를 다루는 데 가장 중요한 공식이지요. 여기서 앞의 er는 복소수의 절대값(크기)이라 하고, theta는 편각이라고 합니다. (복소평면에 그려보면 x축과 이루는 각이 theta, 길이가 er이 나오니까요.)

두 복소수를 곱하면, 복소수의 절대값은 곱해지고, 편각은 더해집니다.


자, 이제 만약에 다음 식을 가지고 생각해볼께요.


아크탄제트 함수를 위처럼 변형해봤습니다. (식1)의 실수부가 분모, 허수부가 분자라는 것을 기억하세요. 다른 말로 쓰면,

는 복소수 가 나타내는 각과 같은 값이라는 것이죠.


이므로, 이 값은 (1 + i)가 나타내는 값과 같습니다.


여기까지 장황하게 설명한 것은 아크탄젠트들의 합차를 복소수의 곱셈으로 표현할수 있음을 설명하기 위해서였습니다. 자, 이제 이 식을 복소수로 바꿔볼까요


->



계산을 해보니 (1+i)의 배수가 나옵니다. 편각을 결정짓는건 실수부와 복소수의 비라는 것을 기억하세요. 편각이 Pi/4와 같다는 알 수 있습니다.


수학자들은 비슷한 기교를 통해 변태..같은 공식들을 많이 만들어냈는데, 다음과 같아요,,,



(핡.. 수학자들은 다들 변태야... 더 끔찍한 식은 여기에서 마음껏 볼수 있습니다.)


모두 다 똑같이 아크탄젠트 급수를 이용한다는 점에서는 변함없지만, 인수가 작아졌기때문에 수렴속도는 훨씬 빨라졌지요.


코드로 구현해봅시다.



먼저 아크탄젠트 함수부터~

저번 연재에서 사용했던 것과 거의 똑같습니다.

이제 이를 이용해 Pi를 계산해보아요.




속도를 한번 비교해볼까요? 길이를 512(512*32비트 정밀도)로 하고, Pi를 구해봤습니다.

pi_machin3함수가 처음구현했던 pi_atan보다 두 배 이상 빠르네요!

(정밀도를 늘리거나 줄여도, 이 비율은 거의 그대로 유지됩니다. 왜냐면 지금까지 구현한 atan, machin1, machin2, machin3 알고리즘은 모두 같은 시간복잡도를 가지고 있거든요.)


하지만 여전히 느립니다. 더 빠른 알고리즘이 없을까요? 다음 연재에서는 훨씬 더 빠른, 코사인 함수를 이용한 incremental한 pi계산법을 다뤄볼게요~

  1. http://en.wikipedia.org/wiki/Machin-like_formula [본문으로]
Tags
이 댓글을 비밀 댓글로
    • 차하나
    • 2013.06.10 01:21
    잘 봤습니다 여쭤볼게있는데요. 머신이 만든 식에서 4arctan(1/5)-arctan(1/239)에서 복소수 -2^2*13^4*(1+i)에 대응되잖아요. 그러면 편각이 5pi/4아닌가요? 모듈로 pi 라서 pi/4인가요? 제가 이것 요즘 엄청 고민해서요. 제가 읽어본 논문에서 arctan(b/a)=arg(a+bi) (mod 2pi)라는데 2pi라는게 도무지 이해가 안됩니다. 예를 들어 a=1, b=1일떄와 a=-1, b=-1일때는 좌변은 같지만 우변은 다르지않나요? 모듈로 2pi가 맞나요? 죄송하지만 마지막으로 한개만 더 여쭤볼게요ㅠㅠ -78125=(278+29i)(2-i)^7 인데요. 그러면 우변은 arctan(29/278)+7arctan(-1/2)이잖아요. 저 함수값 생각해보면 대충 음수나올거라는거 알잖아요. 그래서 함수값은 -pi가 되겠는데, 여기서 고민이되는게 -78125의 편각은 pi도 되고 -pi도 되는데.. 둘 중에 하나 딱 고를수있는 방법이있나요? 얘기가 길어졌네요.. 죄송합니다..
    • 편각을 모듈로 2Pi로 제한하는것은 삼각함수의 주기성 때문에 어차피 2Pi마다 같은값이 돌아오기 때문이죠. 그래서 정의상 편각의 범위는 (-Pi, Pi]가 됩니다. (물론 다르게 잡을 수도 있고, [0, 2Pi)로 잡는 경우도 있고합니다 이건 쓰는 사람마음대로니까요)
      그래서 마지막 질문에 대한 답은 편각 범위를 어떻게 정의했느냐에 맞춰서 고르시면 될거같습니다.
    • arg(a+bi) = arctan(b/a) 라는 편각 계산방법은 a와 b가 모두 양수일때만 사용가능한걸로 알고있습니다. 말씀하신 것처럼 b/a값이 같아도 편각이 다른 경우가 있죠. 그래서 편각을 계산할때는 a+bi가 어떤 사분면에 위치하는지를 보고 계산해야하는 걸로 알고 있어요.
      http://en.wikipedia.org/wiki/Argument_%28complex_analysis%29#Computation
      에 편각을 계산하는 정확한 방법이 나와있으시 참고하시면 도움이 될거 같습니다.