이전글 : 유클리드 호제법(http://www.tibyte.kr/224)



이진 최대공약수 알고리즘(binary GCD algorithm은 스테인 알고리즘(Stein's algorithm)이라고도 하며 비트연산을 이용하여 컴퓨터에서 최대공약수를 (나머지연산을 사용하는 유클리드 알고리즘보다) 빠르게 구할 수 있는 알고리즘이다.


이 알고리즘으로 36과 24의 최대공약수를 구하는 예시를 보면 다음과 같다.

36 → 18 →9 →9 →6 →3 →3 

24 → 12 →6 →3 →3 →3 →0


여기서 이런 규칙을 볼 수 있다.

- 두 수가 짝수일 때는 2로 나눈다.

- 한 수만 짝수일 때는 그 수를 2로 나눈다.

- 두 수가 홀수일 때는 큰 수에서 작은 수를 뺀다.

- 0이 되는 수가 있으면 종료한다.


여기서 최대공약수를 얻는 방법은

두 수가 짝수여서 2로 나누는 계산을 할 때마다, 변수 a를 1씩 증가시킨다.

그리고 마지막에 남은 0이 아닌 수를 b라고 하면

를 계산하는 것이다.


두 수가 모두 2로 a번 나눠진다는 것은

두 수의 공약수에 2가 a번 포함되어있기 때문이라고 생각하면 이해하기가 쉽다.


36과 24의 경우에서, 마지막에 3이 남았다는 것은

바로 전 단계에 두 수에 모두 3이 포함되어있는 것이기 때문에 이 역시 곱해준다(위 식에서는 b)


규칙에서 두 수가 홀수일 때 값을 빼는 것은 유클리드 알고리즘과 관련있는 부분이다.




C언어 코드

int calcGcd(int n1, int n2)
{
    int exp = 0;
    int n;

    while (((n1|n2)&1) == 0) {
        ++exp;
        n1 >>= 1;
        n2 >>= 1;
    }

    while (n1!=0 && n2!=0) {
        while ((n1&1) == 0)
            n1 >>= 1;
        while ((n2&1) == 0)
            n2 >>= 1;

        if (n1 > n2)
            n1 = n1-n2;
        else
            n2 = n2-n1;
    }
    if (n1!=0)
        n = n1;
    else
        n = n2;

    return n<<exp;

}





관련글 : 유클리드 호제법(http://www.tibyte.kr/224)









유클리드 호제법(Euclidean algorithm)으로 최대공약수를 구할 수 있다.


두 정수 m과 n (m>n)의 최대공약수를 구하는 문제를

m-n과 n의 최대공약수를 구하는 문제로 줄일 수 있다.


왜냐하면 m과 n의 최대공약수를 G라 하면 다음과 같이 나타낼 수 있는데

m = Ga

n = Gb

m-n = G(a-b)와 같이 묶을 수 있기 때문이다. (a와 b는 서로소(relatively prime))



예를 들어 94와 30이 있을 경우 아래와 같은 단계를 거쳐 문제를 줄여가면서

최대공약수인 2를 구할 수 있는 것이다.


94    30   (A)

64    30   (A)

34    30   (A)

 4     30   (A) (B)

 4     26   (B)

 4     22   (B)   

 4     18   (B)

 4     14   (B)

 4     10   (B)

 4      6    (B)

 4      2    (B)

 2      2

 2      0


위 단계에서 (A)부분과 (B)부분을  잘 보면 뺄셈 연산을 여러 번 수행한 결과가 결국 나머지 값이라는 사실을 알 수 있다. 따라서 코드로 작성할 때 나머지연산(%)을 사용하면 다음과 같이 단계를 줄일 수 있다.


94   30

 4    30

 4     2

 2     2

 2     0



아래는 C언어로 작성한 코드이다.


int gcd(int n1, int n2)
{
    int temp;
    if (n1 < n2) {
        temp = n1;
        n1 = n2;
        n2 = temp;
    }
    while (n2 != 0) {
        temp = n1%n2;
        n1 = n2;
        n2 = temp;
    }
    return n1;
}






관련글 : 이진 최대공약수 알고리즘 (http://www.tibyte.kr/225)














http://tibyte.kr/218 에서 이어집니다.




주의) Matlab에서는 특정 분포를 하는 다양한 난수생성을 지원합니다.

이 포스트는 코드를 간단히 구현해 보기 위한 목적으로 matlab을 사용합니다.



정규난수를(normal random variable)얻기 위해 박스-뮬러 변환(Box-Muller Transform) 을 쓸 수 있지만


연산과정에 삼각함수가 들어 있기 때문에

이 복잡성을 해결하기 위해 George Marsaglia가 고안한 극좌표법(polar method)을 쓸 수 있다.




박스뮬러 방법에서 균일랜덤변수 u를 사용하여  sin(2pi*u)와 cos(2pi*u)로 나타냈던 것을

극좌표법에서는 원 위에 무작위로 점(u1, u2)을 찍어서 cos이 밑변/빗변,

sin이 높이/빗변인 것을 이용하여  삼각함수를 대체한다. 


박스뮬러 방법에서의 U1이 거리제곱(sq)으로 대체되고,

삼각함수 대신 제곱과 제곱근연산이 들어가게 된다.

[그림 1]





매틀랩 코드로 구현한 결과이다. 표준편차가 1이고 평균이 0인 정규분포를 생성. 


N = 1000000;

u1 = zeros(1,N);

u2 = zeros(1,N);

wsq = zeros(1,N);

for i=1:N

   sq = 2;

   while sq>1

      r1 = rand;

      r2 = rand;

      sq = r1.*r1 + r2.*r2;

   end

   u1(i) = r1;

   u2(i) = r2;

   wsq(i) = sq;

end

s = 1;

m = 0;

x = m+s*sqrt(-2*log(wsq)).*u1./(wsq.^0.5);

y = -1.*(m+s*sqrt(-2*log(wsq)).*u2./(wsq.^0.5));

n = horzcat(x,y);

hist(n, 50);



히스토그램 결과.


[그림 2]





이 방법의 단점은 [그림 1]과 코드에서 볼 수 있듯이

점이 원 밖에 찍혔을 경우 평균적으로 1.27배 만큼 난수를 다시 뽑아야 한다는 것이다.

(1.27배는 기하분포의 급수로 계산)










매틀랩에서는 normrnd()함수로 정규난수(normal random variable)를 발생시킬 수 있긴 하지만

여기서는 균일분포(균등분포)를 변환하여 생성해보도록 한다. 



균일분포(uniform distribution)인 랜덤변수와

박스뮬러 변환(Box-Muller Transformation)을 이용하여 정규난수(Gaussian Random)을 만들 수 있다.


아래는 매트랩으로 작성한 코드와 실행결과이다.



format long

N = 1000000;

u1 = rand(1,N);

u2 = rand(1,N);

s = 1;

m = 0;

x = m+s*sqrt(-2*log(u1)).*cos(2*pi*u2);

y = m+s*sqrt(-2*log(u1)).*sin(2*pi*u2);

n = horzcat(x,y);

hist(n, 50);









- 매틀랩의 다양한 plot들을 사용해서 그려본 도표들.










- 누적분포함수(cdf)






극좌표 방법(polar method)에 대해서는 다음 글에 이어집니다.


관련글 : http://tibyte.kr/219






세 점(기준점, 점1, 점2)이 있을 때 atan()함수를 1번만 써서 사잇각을 구하는 방법이다.

원리는 간단하다.


아래와 같이 한 각과 기준각이 이루는 각도를 α, 나머지 각을 β라 한다.



그러면 θ = α - β 라 할 수 있다.


tanθ는 다음과 같이 나타낼 수 있으므로 sinθ와 cosθ를 구하여 atan2의 각 매개변수로 넣으면 각을 알 수 있을 것이다.


먼저 sinθ는 삼각함수의 덧셈정리에 의해 이렇게 전개할 수 있다.


위 식을 처음에 주어진 x분과 y성분으로 나타낼 수 있으므로,

최종적으로는 이렇게 쓸 수 있다.




코사인세타도 같은 방법으로 전개한다.


공통항인 1/l1l2를 소거하면 아래와 같은 표현식을 얻을 수 있다.



angle = atan2(y1x2-x1y2, x1x2+y1y2);




atan2 함수가 아닌 atan함수를 사용한다면 atan(sin값/cos값) 형태로 사용하면 될 것이다.

→ angle = atan(y1x2-x1y2 / x1x2+y1y2);



▼ 컴퓨터 프로그램으로 작성하여 구동한 화면.






( tibyte.kr과 cafe.naver.com/sdbx에 게시)


 2D 게임을 짜다가 날아가는 미사일이 적을 향해 회전해야 하는 것을 구현할 때

미사일이 시계방향(CW : ClockWise)으로 회전할 것인지, 아니면

반시계방향(CCW : CounterClockWise)으로 회전해야 하는지 알아내야 한다.


 arctangent 함수로 미사일 진행방향과 적(목표물)의 각도를 구하여

미사일의 현재 진행각도와 조건문으로 비교하면 쉽게 될 것 같지만

각도체계가 -pi ~ pi 로 pi각 부분에서 불연속적인 수치가 나타나기 때문에

조건문이 생각보다 복잡해진다.



 미사일 진행방향의 방향벡터를 u, 미사일로부터 목표물에 이르는 방향벡터를 v로 놓고

두 방향벡터의 판별식을 보면 시계방향으로 돌아야 하는지 반시계방향으로 돌아야 하는지 알 수 있다.

판별식 값의 크기는 기하학적으로 두 방향벡터가 이루는 평행사변형의 넓이가 되는데

여기서 그 값의 부호를 보면 두 벡터의 위치관계가 나온다.



[그림1]



[그림1] 과 같이 방향벡터 u를 기준으로, det[u v]<0이면 v는 u보다 시계방향으로 회전한 곳에 위치하고,

det[u v]>0이면 v는 u에서 반시계방향으로 회전한 곳에 위치한다.

판별식 값이 0이면 두 방향벡터는 평행하다.






[그림2]



 다시 처음의 미사일 문제로 돌아와서 [그림2]에는

미사일의 벡터(벡터 AB)와 목표물(점 P)가 그려져 있다.

여기서 벡터AB와 벡터AP를 구하여 판별식을 구해본다면

미사일이 어디로 돌아야 하는지 계산할 수 있을 것이다.


벡터 AB : (Bx-Ax, By-Ay)

벡터 AP : (Px-Ax, Py-Ay) 




det[AB AP] = (Bx-Ax)*(Py-Ay) - (By-Ay)*(Px-Ax)



즉  (Bx-Ax)*(Py-Ay) (By-Ay)*(Px-Ax)의 부호를 보면

시계방향, 반시계방향을 결정할 수 있다는 것이다. 


컴퓨터 그래픽에서는 y좌표의 값이 아래로 갈수록 커지기 때문에

위의 내용을 반전하여,

판별식값이 0보다 작으면 반시계방향이고,

판별식값이 0보다 크면 시계방향으로 처리하면 된다.









1페이지부터 5페이지까지 5개의 장면이 있고

다음 버튼과 이전 버튼이 있어서 페이지를 넘길 수 있는 프로그램을 만드는데

if문을 쓰지 않고 모듈러연산으로 처리하는 것이 생각나서 식을 구해보기로 했다.


아래 그림은 1부터 n까지의 자연수가 있을때

아무것도 더하지 않고 mod n연산을 하여 +1을 한 결과와,

처음에 1을 더한 경우, 2를 더한 경우에 대한 내용이다.


아래 그림은 처음에 n-2을 더한 경우와, n-1을 더한 경우이다. 

예를 들어 자연수의 집합 {1, 2, 3, 4, 5} 가 있을 때

집합의 각 원소에 n-2인 3을 더하면  집합은{4, 5, 6, 7, 8}가 되고 

이 집합의 mod 5를 구하면 {4, 0, 1, 2, 3} 이 되어

다시 이 집합에 +1을 하면 {5, 1, 2, 3, 4}가 되어서

최종적으로는 원래 집합이 우측으로 1칸씩 순환 이동하였다고 볼 수 있다.

또한 2는 1이 되고, 3은 2가 되고, 4는 3이 되고, 5는 4가 되고, 1은 다시 5가 되었으므로

처음에 만들고자 하였던 프로그램의 이전 버튼에 사용할 수 있을 것이다.


다시 위의 두 그림을 자세히 살펴보면,

집합의 원소 a가 있고 집합의 원소가 n개인 경우에, 이 집합의 원소를 p씩 증가시키고 싶을 때

(예를 들어 각각의 1,2,3,4,5 페이지에서 다음버튼을 눌러서 각각 2,3,4,5,1 페이지를 만들고자 하는 경우 등) 

다음과 같은 함수를 얻을 수 있다.

f(a,p,n) = (a+(p-1))% n +1

여기서 p는 좌측으로 p칸 순환시프트시킬 값도 된다.


처음의 프로그램 문제로 되돌아가서

5페이지일때의 다음버튼을 만든다고 하면 아래와 같이 프로그램을 작성할 수 있을 것이다.

page = (currentPage+(1-1))%5 + 1;

필요없는 연산을 제거하면,

page = currentPage%5 + 1;


이번에는 2페이지씩 넘기는 버튼을 만든다고 하면 코드는,

page = (currentPage+(2-1))%5 + 1; 에서 역시 필요없는 부분을 제거하여
page = (currentPage+1)%5 + 1;
이 된다.

그러면 반대방향으로 1페이지를 넘기는 이전버튼에 들어갈 코드는 어떻게 될까?
순환 시프트이므로 p값을 넣을 자리에 n-p값을 놓으면 된다.
page = (currentPage+(5-1-1))%5 + 1; 
page = (currentPage+3)%5 + 1; 


정리 :
(a+(p-1))% n +1
→ p씩 순환증가
→ p씩 좌측 순환시프트 

(a+(n-p-1))% n +1
→ p씩 순환감소
→ p씩 우측 순환시프트 



================================================================================



+ Recent posts