나의 큰 O는 log x야

저번 연재에서 아크탄젠트 급수를 이용해서 Pi를 계산해보았는데, 머리를 조금 더 굴리면 훨씬 더 빠른 알고리즘을 만들어낼 수 있습니다. 가장 대표적인 공식은 다음과 같아요. John Machin이 1706년에 발견한 급수인데, Machin의 이름을 따서 Pi를 계산해내는 일련의 아크탄젠트 공식들을 Machin-like 공식이라고 부릅니다. 그냥 지나가면 재미없으니, 저 공식이 어떻게 유도되었는지 거슬러 올라가봅시다. (복소수를 이용한 증명을 소개하니, 이해하려면 복소수에 대한 기초적인 지식이 필요해요..) ...(식1) 복소수 지수를 다루는 데 가장 중요한 공식이지요. 여기서 앞의 er는 복소수의 절대값(크기)이라 하고, theta는 편각이라고 합니다. (복소평면에 그려보면 x축과 이루는 각이 theta..

프로그래밍/Multi precision 2012. 11. 22. 01:38

드디어! 이번 연재의 궁극적인 목표인 Pi 계산에 필요한 모든 것들이 준비되었군요. (사실 Pi를 계산하기 위해서는 사칙연산만 구현하면 충분합니다...ㅋ) 미적분이나 급수가 등장하기 전까지 고대인들은 Pi를 구하기 위해서 여러가지 삽질...(보다는 노력!)을 했었는데요, 수학적 테크닉이 발전하면서 빠른 속도로 Pi를 구하는 방법이 어마어마하게 등장했습니다. 이제 우리는 책상머리에 앉아서 몇개의 수식을 가지고 Pi를 계산해 낼 수 있지요. 가장 널리 쓰이는 방법은 Machin의 공식입니다. 이 공식을 본 수많은 사람들은, Pi가 이렇게 단순하게 홀수들의 역수의 합차로써 계산된다는것에 충공깽했죠. 이 식은 아크탄젠트 급수에서 쉽게 유도될수 있습니다. (삼각함수 정도는 기본인거잖아요? 삼각함수/라디안에 대한 ..

프로그래밍/Multi precision 2012. 11. 17. 15:45

최적화하는건 머리가 너무 아프므로, 이번엔 산뜻하게 제곱근 계산하는 함수나 만들어보렵니다. 사칙연산을 다 구현하면 BigFloat 자료형을 다 만들었다고 착각하면 아니되옵니다.ㅋ 이제 시작이지요. 그중에서 제일 기본적인게 제곱근입니다. 계산도 어렵지 않아요. 나눗셈 구현할때 사용했던 Newton-Raphson 방법을 기억하시려는지요? 이를 제곱근 계산에도 이용할 수 있습니다.(바빌론에서 제일 처음 사용되었다고 해서 바빌로니아 방법이라고 부르기도 합니다.) 이 공식 기억나지요? f(x)=0의 근을 근사해주는 Newton-Raphson의 공식이었죠. 우리는 제곱근을 구하고 싶으므로 f(x)를 다음과 같이 놓겠습니다. (f(x)의 근은 sqrt(S) (혹은 -sqrt(S) )가 되겠지요?) 그럼 아래와 같이 ..

프로그래밍/Multi precision 2012. 11. 8. 21:28

요즘 컴파일러들은 매우 똑똑합니다. 어지간한 코드는 알아서 최적화시켜주지요. 그래서 대부분의 경우에는, 힘들여 어셈블리어 코드를 써서 최적화해놔도 컴파일러에 의한 최적화보다 못한 경우가 많습니다.하지만! 이렇게 큰 수끼리 연산을 하는 코드는 최적화할 가능성이 많습니다. 우리가 예전에 덧셈과 뺄셈을 어떻게 구현했는지 한번 떠올려볼까요?(2. 덧셈, 뺄셈 구현하기) 긴 정수의 덧셈을 다음과 같이 구현했었죠.1. 먼저 두 정수를 더한다.2. 더한 결과가 원래 정수보다 작으면 받아올림이 발생한거다.3. 윗 자리를 더할때 받아올림된 1도 같이 더해준다. 이를 위해서 n자리의 정수라면 총 n번의 덧셈과 n번의 비교와, 받아올림이 있다면 n번의 추가적인 덧셈이 필요했습니다. 즉 최악의 경우 3n번의 덧셈/비교 연산이..

프로그래밍/Multi precision 2012. 11. 5. 00:54

저번 연재에서 나눗셈 구현한다고 그러면서 계속 역수 계산하는것만 구현하고 정작 나눗셈은 구현하지 않았네요ㅋ. 간단하게 보충하고자 합니다.역수 계산하는 함수를 구현한 이상 나눗셈을 구현하는것은 어렵지 않습니다. 그냥 역수를 곱하면 되는것이니깐요.코드를 볼까요? int RR_div(PRR dest, PRR a, PRR b) { if(RR_iszero(b)) return -1; PRR rcp = RR_create(dest->size 1); RR_inv(rcp, b); RR_mul(dest, a, rcp); RR_delete(rcp); return 0; } 참... 쉽죠?역수도 구하고, 곱셈도 해야하므로, 당연히 곱셈연산보다는 나눗셈 연산이 조금 느릴수밖에 없습니다..ㅠ 하지만 시간복잡도는 곱셈과 마찬가지로 O..

프로그래밍/Multi precision 2012. 11. 4. 23:49

전 연재에서 Newton-Raphson 방법으로 역수를 계산해보았는데요, Newton-Raphson 방법의 특징을 이용해서 코드를 조금더 최적화할 수 있습니다. Newton-Raphson 방법은 수치 계산에 있어서 점진적(incremental)입니다. 무슨 말인가 하면, 반복적으로 같은 알고리즘을 적용하여 나가면서 수치 계산의 정밀도를 높여가는데, 중간에 계산을 멈추었다가도 그 중간결과를 이용해 다시 계산을 시작해 나가면 더 높은 정밀도를 얻을 수 있다는 것입니다. Goldschmidt의 방법은 연속적으로을 곱해나가면서 역수를 계산해내었는데요, 이는 의 값을 따로 알고있어야만 다음 과정으로 넘어가서, 더 높은 정밀도의 결과를 얻을수 있지요. 즉 이 알고리즘은 incremental하지 않습니다. incre..

프로그래밍/Multi precision 2012. 11. 1. 14:19

저번 연재에서는 Goldschmidt의 방법을 이용해서 나눗셈을 구현했었는데요, 이번에는 Newton-Raphson의 방법을 이용해서 구현해보려고 합니다. 그리고 그 둘 중에서 어떤 놈이 빠른지 확인해보는거지요. Newton-Raphson 방법에서는 초기값이 중요합니다. 초기값을 잘못 잡으면 아예 수렴 방향이 바뀌어서 엉뚱한 곳으로 수렴할 수가 있거든요. (기억이 나지 않는다면 BigFloat로 Pi를 구해보자-4.나눗셈은 어려워 (나눗셈 구현 알고리즘) 참고) 초기값을 적당히 잡아서 i를 키워나가면 Xi는 1/D에 수렴하겠지요. 그러면 초기값은 어떻게 잡을까요? D의 범위를 0.5 < D < 1로 잡았으니, 1/D의 범위는 1 < 1/D < 2 가 될겁니다. 그냥 단순하게 1차 선형으로 근사를 시켜보면..

프로그래밍/Multi precision 2012. 10. 28. 22:57

저번 연재에서 나눗셈 알고리즘을 설명했으니, 이번에는 곧바로 구현에 들어갑니다. 먼저 나눗셈을 위해서는 비트 시프트 연산이 필수적입니다. 다음은 16비트 정수를 왼쪽으로 3비트 시프트하는 연산입니다. a0, a1, a2는 오버플로우로 버려지게 되고, 오른쪽의 비는 3자리는 0으로 채워지게 됩니다. 다음은 오른쪽으로 3비트 시프트하는 연산입니다. b5, b6, b7은 언더플로우로 버려지게 되고, 왼쪽의 비는 3자리는 0으로 채워지게 됩니다. 자 그럼 이제 코드를 볼까요? /* len길이의 dest 배열을 왼쪽으로 shift 비트만큼 시프트하는 함수. */ inline void i_shls(RRU32* dest, RRU32 shift, RRU32 len) { RRU32 w = shift >> 5; shift..

프로그래밍/Multi precision 2012. 10. 24. 23:06