이전글 : 유클리드 호제법(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보다 크면 시계방향으로 처리하면 된다.







 







점(a,b) 와 점(c,d)가 주어졌을 때, 정삼각형이 되기 위한 나머지 한 점은
P1과 P2. 2개가 존재합니다. (a,b 점과 c,d점 사이의 거리를 반지름으로 하는 두 원을 그리면 그 두 원의
교점이 P1과 P2가 됩니다.)

위의 그림과 같이 중점을 m이라 놓고 P1 = (x1,y1) , P2 = (x2,y2)라고 놓습니다.
 

그리고 아래와 같이 직각삼각형을 그렸을 때 선분PM 과 주어진 두 점을 잇는 선분이 수직이므로
아래 두 삼각형은 닮음이 됩니다.



아래쪽 삼각형에 대한 위쪽 삼각형의 닮음비가
2분의 루트3 대 1 이므로,

P1점의 x좌표는, m의 x좌표에다가 (d-b)*닮음비가 되고,
y좌표 역시 m의 y좌표에 (c-a)*닮음비가 됩니다.




P2점의 좌표도 같은방법으로 구합니다.



아래와 같이 정리가 됩니다.






또 다른 방법은
역시 두 점(a,b),(c,d)가 주어졌을때,
 다각형(x1,y1)(a,b)(c,d)가 정삼각형이므로 사잇각이 60도가 됩니다. 



여기에 삼각함수의 덧셈정리로 주어진 선분의 각도 세타를 소거해주면
좌표를 하나씩 구해나갈수 있습니다.








copyright(C) www.tibyte.kr

+ Recent posts