고정 헤더 영역
상세 컨텐츠
본문
저번 연재까지는 아크탄젠트 급수와, 이를 이용한 Machin 급수로 Pi를 계산하는 방법에 대해서 다뤘습니다. 이번에는 저번 방법보다 훨씬 빠른 코사인 함수를 이용한 계산법에 대해서 소개합니다.
위 점화식은 n이 무한히 커짐에 따라 Pi/2에 수렴합니다. x 값의 정밀도가 어느 정도인지에 상관없이 항상 그보다 더 정밀한 근사값을 계산해줍니다.
오차가 얼마나 빨리 줄어드는지 계산해 볼까요?
오차를 입실론이라고 하면 다음과 같이 식을 세울 수 있습니다.
함수 꼴로 바꿔보면 조금더 보기 쉽겠네요
테일러 급수를 이용하면 다음 사실이 자명함을 알 수 있습니다.
그러므로 역시 다음도 성립합니다.
즉 오차는 매 반복마다 세제곱 한것보다 더 줄어듭니다. 다르게 말하면, 매 반복마다 정밀도가 최소 3배씩은 증가한다는 말입니다.
수렴속도도 빠를 뿐만 아니라, incremental 하다는 장점도 있습니다. (나눗셈 구현할때 사용했던 Newton-Raphson 방법처럼ㅋ : 7. 나눗셈 최적화와 큰O 표기법(Big-O Notation))
즉 처음에는 작은 정밀도로 계산해나가다가, 나중에 차차 정밀도를 키워가는 테크닉을 사용할 수 있다는 것이지요.
이 방법을 구현하기 위해서는 먼저 코사인 함수를 구현해야 합니다. 코사인 함수의 테일러 전개는 널리 잘 알려져 있지요.
분모가 팩토리얼로 늘어나면서 x값이 무엇이든간에 상관하지 않고 다 수렴시켜버립니다!! (오오미 팩토리얼 증가의 위대함!)
하지만 이 전개가 모든 값에 다 수렴한다고 해도, 빠른 수렴속도를 위해서는 최대한 작은 x를 고를 필요가 있습니다. 그래야 수렴속도가 빨라지니까요. 우리가 x로 사용할 값은 Pi/2 근처값이니 1보다 큽니다. 1보다 작으면 수렴 속도가 빨라서 참 좋을텐데요...
일단 코사인 함수를 구현해봅시다!
매번 x^2을 계산하는것은 시간/에너지 낭비니깐, 먼저 x^2을 계산해서 저장해두면 좋겠네요.
그리고 처음 계산한 x^2/2!는 다음에 x^4/4!를 계산할때 이용하면 좋겠네요. x^2/2! 에다가 x^2을 곱하고, 3*4를 나누면 바로 다음 항을 구할 수 있으니까요. 이런 테크닉을 이용해 3번의 곱셈과 1번의 덧셈을 이용해 하나의 항을 계산할 수 있습니다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 |
void RR_cos(PRR dest, PRR a){//dest = cos(a) PRR a2 = RR_create(dest->size + 1); PRR t = RR_create(dest->size + 1); PRR u = RR_create(dest->size + 1); PRR rcp = RR_create(dest->size + 1); RR_initu(dest, 1, 0); //dest = 1 RR_initu(t, 1, 0); //t = 1 RR_mul(a2, a, a); //a2 = a^2 for(RRU32 i = 1; i < 0xFFFFFFFE; i+=2) { RR_invu32(rcp, i*(i+1)); RR_mul(u, t, a2); RR_mul(t, u, rcp);//t = t* a^2 / (i*(i+1)) if(RR_iszero(t))break;//t가 0이 되면 탈출 RR_neg(t); RR_add(dest, t); } RR_delete(rcp); RR_delete(t); RR_delete(u); RR_delete(a2);} |
불행히도 우리는 저 반복문이 몇번이나 돌지를 사전에 예측할 수가 없습니다 ㅠ(O(n)보다는 작지만, O(log n)보다는 큰 그 사이 어딘가에...). 확실한 건 초기값 a가 작아서, 항이 빠르게 감소하길 기대하는 것 뿐이지요. 곱셈이 시간복잡도가 O(n2)이었으니 코사인함수 전체의 시간복잡도는 O(n2log n)보다는 크고 O(n3)보다는 작게 되겠습니다.
하지만 저 속도에는 만족못하겠네요. 삼각함수의 반각공식이 있습니다.
cos x 를 계산하는 대신에, cos x/2를 계산해서 이것으로 cos x를 구하는 겁니다. 이 방법을 재귀적으로 사용하면 빠른 속도로 cos 을 계산해 낼수 있을 겁니다. x를 너무 많이 반으로 쪼개어 계산하면 오차가 누적되어 다른 값이 나올 수 있으니 주의해야합니다.
이 방법을 이용한 cos2함수입니다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 |
void RR_cos2(PRR dest, PRR a){ PRR ah = RR_create(dest->size + 1); PRR dh = RR_create(dest->size + 2); PRR eh = RR_create(dest->size + 2);//오차 방지를 위해 넉넉히 공간을 잡는중 RR_shrcopy(ah, a, 32);//2로 32번 나눔 RR_cos(dh, ah);//위에서 구현한 cos함수 호출 for(int i = 0; i < 31; ++i) { RR_mul(eh, dh, dh);//제곱하고 RR_shl(eh, 1);//2를 곱하고 if(i_subu32(eh->frac, eh->size, 1, 0))//1을 빼준다 { eh->sign = ! eh->sign; i_cmplt(eh->frac, eh->size);//부호가 바뀌었을 경우 처리 } PRR t = eh; eh = dh; dh = t; } RR_mul(dest, dh, dh); RR_shl(dest, 1); if(i_subu32(dest->frac, dest->size, 1, 0)) { dest->sign = ! dest->sign; i_cmplt(dest->frac, dest->size); } RR_delete(eh); RR_delete(dh); RR_delete(ah);} |
RR_cos2함수는 RR_cos함수보다 훨씬 빠릅니다. 문제는 약간 오차가 발생한다는 건데, 그 정도는 충분히 감수할 만하지요. 오차문제는 dest의 size를 처음에 넉넉하게 잡음으로써 해결할 수 있습니다.
자 이제 우리의 최종 목표였던 Pi를 계산해봅시다.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 |
void RR_pi_cos(PRR dest){//RR_cos 함수를 사용하는 버전 PRR t = RR_create(dest->size); RR_init(dest, 1.57); t->size = 3; dest->size = 3; RRU32 c = 123456789; for(int i = 0;; ++i) { RR_cos(t, dest); RR_add(dest, t);//dest += cos(dest) if(t->size == t->alloced) { if(c == dest->frac[dest->size - 2])break;//수치가 안정화되면 탈출 c = dest->frac[dest->size - 2]; } t->size = MIN(t->size*3, t->alloced);//정밀도를 3배씩 증가시킵니다. dest->size = t->size; } RR_shl(dest, 1);//결과물 2배해주기 RR_delete(t);}void RR_pi_cos2(PRR dest){//RR_cos2를 사용하는 버전 PRR t = RR_create(dest->size); RR_init(dest, 1.57); t->size = 3; dest->size = 3; RRU32 c = 123456789; for(int i = 0;; ++i) { RR_cos2(t, dest); RR_add(dest, t);//dest += cos(dest) if(t->size == t->alloced) { if(c == dest->frac[dest->size - 2])break;//수치가 안정화되면 탈출 c = dest->frac[dest->size - 2]; } t->size = MIN(t->size*3, t->alloced);//정밀도를 3배씩 증가시킵니다. dest->size = t->size; } RR_shl(dest, 1);//결과물 2배해주기 RR_delete(t);} |
32비트 512 길이 (즉 32*512비트 정밀도)
RR_pi_cos2 함수의 압도적인 승리입니다!
'프로그래밍 > Multi precision' 카테고리의 다른 글
| [진짜 BigFloat 구현하기] 0. 시작에 앞서서 (0) | 2013.11.20 |
|---|---|
| 환형 큐를 사용한 다중정밀도 연산 (0) | 2013.10.26 |
| BigFloat로 Pi를 구해보자-13. 가우스-르장드르 알고리즘 (6) | 2012.12.01 |
| BigFloat로 Pi를 구해보자-11. Pi 구하기: Machin 급수 (3) | 2012.11.22 |
| BigFloat로 Pi를 구해보자-10. 본격 Pi 구하기: 아크탄젠트 급수 (0) | 2012.11.17 |
| BigFloat로 Pi를 구해보자-9. 제곱근 구현하기 (0) | 2012.11.08 |
댓글 영역