레이블이 수학인 게시물을 표시합니다. 모든 게시물 표시
레이블이 수학인 게시물을 표시합니다. 모든 게시물 표시

2009년 9월 6일 일요일

에라토스테네스의 체가 과연 빠르긴 빠르네

꼭 이런 짓을 하고싶을 때가 있다.
소수의 합을 구할 때 에라토스테네스의 체가 빠르다는 거 당연한데, 굳이 일일이 계산하는 거랑 비교해보고 싶었다.
왜 그런지 따윈 없고... 단지 있다면 얼마 전 모 블로그에 내가 쓴 답글이 신경쓰여서랄까나...

그래서 VS .Net 2003으로 만들어봤다.

#include "stdafx.h"
#include <memory.h>
#include <math.h>
#include <windows.h>

#define PRIMES 2000000

int is_prime(int a)
{
    int i, sqrn=(int)sqrt((double)a);

    for(i=3; i<=sqrn; i+=2)
        if(a%i==0)
            return 0;
    return 1;
}

void brutal(int &iPrimes, long long &lSum)
{
    int i;
    lSum=2;
    iPrimes=1;

    for(i=3; i<PRIMES; i+=2)
    {
        if(is_prime(i))
        {
            iPrimes++;
            lSum+=i;
        }
    }
}

void eratosthenes(int &iPrimes, long long &lSum)
{
    bool *bNoPrime = new bool[PRIMES];
    lSum = 0;
    iPrimes = 0;
    int iMaxCompare = (int)sqrt((double)PRIMES);

    memset(bNoPrime, 0, sizeof(bool)*PRIMES);

    int i;
    for (i=2; i<=iMaxCompare; i++)
    {
        if (!bNoPrime[i])
        {
            iPrimes++;
            lSum += i;
        }

        for (int j=i*i; j<PRIMES; j+=i) bNoPrime[j] = true;
    }

    for (i=iMaxCompare+1; i<PRIMES; i++)
    {
        if (!bNoPrime[i])
        {
            iPrimes++;
            lSum += i;
        }
    }

    delete bNoPrime;
}

int main(int argc, char* argv[])
{
    DWORD dw;

    int iPrimes;
    long long lSum;

    dw = GetTickCount();
    eratosthenes(iPrimes, lSum);
    dw = GetTickCount()-dw;

    printf("Sieve of eratosthenes method\n  : %d prime numbers, sum is %I64d (%u milisec)\n", iPrimes, lSum, dw);

    dw = GetTickCount();
    Sleep(100);
    brutal(iPrimes, lSum);
    dw = GetTickCount()-dw;
    printf("Brutal force method\n  : %d prime numbers, sum is %I64d (%u milisec)\n", iPrimes, lSum, dw);

    return 0;
}

결과는 당연히 에라토스테네스의 체 쪽이 훨씬 빠르다.

Sieve of eratosthenes method
  : 148933 prime numbers, sum is 142913828922 (31 milisec)
Brutal force method
  : 148933 prime numbers, sum is 142913828922 (469 milisec)

이걸 해보며 발견한 시궁창같은 현실.

1. Visual C++ 6.0은 long long 형도 인식하지 못한다. 얜 정말 안 되는 애다.
첨엔 결과가 무식하게 클 줄 알고 long long으로 정의했는데, long long을 인식못하는 걸 보고 걍 .Net 2003으로 변절.
__int64로 바꾸면 되긴 하지만, 이미 난 삐졌음!
그런데, 결과는 long int 범위 안쪽. OTL

2. VS 계열에선 bool은 내부적으로 1바이트를 사용하고, BOOL(==int)은 4바이트를 사용한다.
따라서 BOOL이 훨씬 빠른 것이 상식이다.
(32비트 CPU에겐 32비트 연산이 가장 빠름, 8비트는 32비트로 확장해서 처리)
그런데, 소수 테이블을 bool로 지정할 때가 BOOL로 지정할 때보다 훨씬 빠르다. 뭐지?

3. VS6이나 VS.Net 2003이나 2백만 개의 배열은 못 잡는다. 런타임 오류 발생.
그래도 new를 이용해 동적으로 할당받는 건 문제가 없으니 예쁘게 봐주고 넘어감.

4. 댓글 보고 다시 확인해보니 결과가 틀렸음. long이 아니라 long long이 맞음.
그런데 잘못 생각한 이유는 다름아닌 .Net 2003의 printf에서 "%lld"를 인식하지 못하기 때문.
MSDN을 뒤져보니 VS2005부터는 %lld를 인식하나보다. 제길슨.


참, 무식하게 찾는 쪽 코드는 Studying the Logical World에서 거의 그대로 업어왔다.
여기가 앞에서 언급한 그 블로그다.

2009년 7월 25일 토요일

등각항로의 정확한 방위/거리를 구하는 방식 발견

앞에 포스팅했던 지구 타원체에서 두 점간의 방위와 최단거리 계산하기를 쓰면서 프로그램을 하나 만들었다.
지표면 상에서 두 좌표 간의 방위와 거리를 각각 대권항로등각항로로 계산해주는 프로그램이다.

그런데, 문제가... 대권항로는 WGS84 좌표계에서 계산을 하는데, 등각항로는 사실 WGA84 타원체의 평균반지름을 반지름으로 갖는 구체에서 계산한다는 것이었다.
(이렇게 계산하면 최악의 경우 0.55%까지 거리의 오차가 발생한다)

오늘... 드디어... WGS84 타원체 상에서 등각항로의 방위와 거리를 정확하게 구하는 방식을 발견했다.
(아쉽게도 논문 같은 것은 찾지 못했음)

Juhani Kaukoranta님께서 Navigation, Astronavigation, Positioning and Geodesy라는 페이지에 식을 포함한 각종 자료를 올려주셨더라.

역시 수학은 멋지고 위대하다.
이 방식을 공부해야겠다!!!

사용자 삽입 이미지

친절하기 짝이 없는 설명의 일부 ^^;


2008년 1월 19일 토요일

Fractal graphics contest 결과 사이트들

사용자 삽입 이미지

FA 2007 : Royalty

수학은 아름답다 : Mandelbrot set이란?에서 프랙탈에 대해서 간단하게 소개했습니다.

FractintUltra Fractal과 같은 좋은 프로그램을 사용하면 기본적인 프랙탈을 쉽게 볼 수 있습니다.
하지만, 좀 더 수준 높은 프랙탈을 보기를 원하신다면 프랙탈 컨테스트 결과를 참고하면 됩니다.

이름이 알려진 프랙탈 컨테스트는 Fractal Art Contest, Ultra Fractal Contest 등이 있습니다.

각 컨테스트의 연도별 수상작을 바로 보려면 다음의 링크를 바로 클릭하시면 됩니다.


정말 정말 아름답지 않나요?

2008년 1월 17일 목요일

수학은 아름답다 : Mandelbrot set이란?

Fractal graphic(이하 프랙탈)의 집합에는 Mandelbrot, Jilua, Newton, Lorenz 등 수많은 집합이 있습니다.
이런 유명한 집합들은 간단한 방식(짧은 수식 한두 개로써!)으로 아름답고 화려한 그림을 보여줍니다.

프랙탈은 자연계에 존재하는 무한반복의 특성을 수식으로 표현하는 방식입니다.
이러한 무한반복 특성을 순환성(Recursiveness)자기 유사성(Self-similarity)으로 정의합니다.
너무 어렵다구요? 아래의 그림들을 보시면 그렇게 어렵지만은 않다는 것을 알 수 있을 것입니다.


사용자 삽입 이미지

Fractal, himself

Benoît Mandelbrot라는 박사님이 계십니다.
이 분께서 바로 Fractal이라는 용어를 만드신 분이십니다. (박수 짝짝짝)

이 분께서 1975년 어느 날 IBM에서 일하시면서 다음의 간단한 수학식을 놓고 연구하고 계셨습니다.

Z = Z²+C

이 수식이 바로 Mandelbrot Set입니다.

복소평면의 각 좌표 Z에서 위의 수식을 무한히 반복하였을 때 각 좌표의 점이 발산하는가 수렴하는가 그리고, 수렴한다면 얼마로 수렴하는가를 복소평면 위에 그려봤는데, 의외의 결과가 나온 것입니다.
(물론, 이 분께는 의외가 아니었겠죠. 평범한 사람들에게 의외란 뜻입니다)

Mandelbrot의 기본 Set는 실수 -2.5~1.5, 허수 -1.5~1.5의 복소평면 위에 그려집니다.
즉, 아래와 같은 복소좌표계에서 각 점에 대해서 같은 수식을 계속 반복한 결과를 색으로 표시한 것입니다.


y=ax+b가 아닙니다. 복소평면입니다.

Mandelbrot set에 대해서 좀 더 궁금하신 분은 Wikipedia를 읽어보시기 바랍니다.





1988년, 80386 시절에 Stone Soup Group에서 Fractint라는 프로그램을 발표합니다. (Fractint 홈페이지)
(정확히는 FRACT386을 발표하고, 1989년에 Fractint로 개명했습니다)

Stone Soup라는 말의 어원은 동화에 나오는 "돌국"이야기 입니다.
어느 지방에 기근이 들었을 때 지나가던 나그네가 돌을 꺼내 국을 끓이면서 마을사람들이 숨겨놓은 재료들을 하나씩 집어넣도록 설득해서 결국 모든 마을사람들이 맛있는 국을 먹을 수 있었다는 이야기가 이름의 어원입니다.
즉, 사람들의 지혜를 모아서 좋은 프로그램을 만들자는 뜻을 담고 있습니다.

당시만 해도 coprocessor가 많이 장착되지 않았습니다. 이 프로그램은 이 사실을 고려해서 가능한 모든 허수 계산을 정수로 계산하도록 만들어졌습니다.
DOS에서 동작하던 이 프로그램은 당시의 척박한 환경에서 정말 빠르게 동작하는 프랙탈 그래픽 도구로 명성을 날립니다. 또한 오픈소스로 개발되어 여러 환경으로 이식되었습니다. 물론, Windows로도 이식되었습니다.

그리고, 평면에 대해서는 좀 더 빠른 알고리즘을 연구하였고 특히, Boundary pass(단축키 x-b)는 특이한 알고리즘으로 대단히 빠른 결과를 보여줬습니다.

Fractint로 Newton set(degree=3)을 그리는 과정 (boundary pass)


이 프로그램(정확히는 DOS용 Fractint)은 지금 사용하기에도 무리가 없는 프로그램입니다.
명확한 그리고, 사용에 정확하게 최적화된 인터페이스도 일품입니다.
하지만, 약간의 단점이 있습니다.

  1. 화려한 메뉴 기능이 없다 (사용에는 전혀 불편하지 않습니다)

  2. Vista에서 동작하지 않는다 (동작하는 모습을 보려면 DOSBox를 사용하면 됩니다)

  3. 기본 설계가 16/256 color를 기준으로 되어있다.



최근, 20.0 발표 이후로 추가적인 개발 소식이 없어서 아쉬웠는데, 그 사이에 Ultra Fractal이라는 상용 프로그램이 나와있었습니다. 2008년 1월 현재 4.04까지 나와있고, 홈페이지에서 free trial 버전도 다운받을 수 있습니다.

1달의 사용제한기한이 있기는 하지만, Fractint와 비교했을 때 엄청난 성능을 보여줍니다.

Mandelbrot set (default)
Newton set (degree=5)


정확하게 같은 수식을 같은 범위에 대해서 같은 크기(1024x768)로 출력한 결과입니다.
정말 멋있고 아름답지 않나요?

마지막으로, 맨 위에 Fractint의 boundary pass를 설명하기 위해 표시한 Newton set(degree=3)의 Ultra Fractal 버전은 아래와 같습니다.

사용자 삽입 이미지
Ultra Fractal로 그린 Newton set(degree=3)



수학이 지루하고, 재미없으며, 그래프가 보기 싫다면 Fractal graphic으로 눈을 돌려보시기 바랍니다.
단순한 수식 1~2개가 너무나 아름답고 환상적인 그림으로 변하는 것을 볼 수 있습니다.

p.s. 사실, 이런 개념을 최초로 생각한 것은 Mandelbrot가 아니라 Gaston Julia라는 프랑스 수학자입니다.
       이 분께서는 무려 1918년에 유사한 개념을 생각했습니다.
       그 시기엔 (그 분께는) 컴퓨터가 없다는 것이 문제였죠. (Mandelbrot 박사님은 IBM 연구원…)



2007년 11월 16일 금요일

사각형의 무게중심 계산법

삼각형의 무게중심 구하는 법은 다들 아시죠?

각 꼭지점에서 마주보는 변의 중점으로 선을 그엇을 때 만나는 점입니다.
좌표(x1,y1)-(x2,y2)-(x3,y3)으로 표현된 경우에는 ((x1+x2+x3)/3. (y1+y2+y3)/3)으로 간단하게 계산할 수도 있습니다.

특징은 중학교 수학책을 참고하시고...

갑자기 궁금해져서 사각형의 무게중심을 구하는 방법을 찾아봤습니다.
사용자 삽입 이미지


의외로 어렵더군요.
오른쪽 그림과 같이 사각형을 두 개의 삼각형으로 쪼갭니다.
(물론 서로 마주보는 두 꼭지점을 연결해야죠)

그러면 삼각형 둘로 나뉘는데,
각각 삼각형의 무게중심을 구합니다.

왼쪽에 있는 삼각형을 A, 오른쪽을 B라고 하면 오른쪽과 같은 그림이 나오게 됩니다. GA는 A의 무게중심, GB는 B의 무게중심입니다.
다음은 GA-GB를 잇는 선분을 면적의 역수로 내분한 점을 찾으면 됩니다.

즉, A의 면적이 4, B의 면적이 1이라고 하면 GA-GB1:4로 내분한 점이 사각형의 무게중심이 됩니다. (4:1이 아닙니다. 그렇게 되면 큰 삼각형 쪽으로 중심이 잔뜩 쏠릴 겁니다)

역시 수학의 세계는 오묘하기 짝이 없습니다~

2차원에서 임의의 한 점이 삼각형 내부에 있는지 판별하는 방법

프로그램 하나를 쪼물딱러리며 만들고 있는데,
임의의 한 점이 삼각형 내부에 있는지 판별하는 방법이 궁금해졌습니다.

이리저리 복잡하게 판단할 수도 있겠지만, 수학적으로 아름다운 방법이 없을까 찾아봤습니다.

그러다가 http://www.mathlove.org/pds/mathqa/faq/geometry/geometry49.html 에서 아래와 같은 답을 찾았습니다.

A=(x_1, y_1), B=(x_2, y_2), C=(x_3, y_3)라고 하고 주어진 점을 (a, b)라고 합시다. 그러면 삼각형의 내부는 아래 세 영역의 교집합입니다.

1. 직선 AB에 의해서 나뉘어 지는 두 영역 중에서 점 C가 속한 영역
2. 직선 BC에 의해서 나뉘어 지는 두 영역 중에서 점 A가 속한 영역
3. 직선 CA에 의해서 나뉘어 지는 두 영역 중에서 점 B가 속한 영역

예를 들어 1의 조건을 식으로 써보면, 직선 AB의 식은

f(x,y) = (x-x_1)(y_1-y_2)-(y-y_1)(x_1-x_2) = 0

이므로 점 C와 P가 같은 영역에 속하는지의 여부는 값

f(x_3, y_3)f(a, b)

의 부호를 살펴보면 됩니다. 양수이면 같은 영역, 음수이면 서로 다른 영역입니다. 이런 식으로 세번의 부등식 판별을 거치면 됩니다.


그 밑에는 간단하게 작성된 C 코드도 있더군요. 대단대단~
그런데, 조금 더 찾아보니 http://gpgstudy.com/forum/viewtopic.php?t=15797&sid=b280649b74e383b21261d6ec8426c5e1 에서 아래와 같은 어마어마한 내용을 찾았습니다.

삼각형이 v1, v2, v3 세점으로 이루어져 있을때
d1 = v2 - v1
d2 = v3 - v1

테스트 하고자 하는 점을 k 라고 할때
p = k - v1

그려면 2차원 상의 한 점 p는 d1과 d2의 선형 합으로 나타낼 수 있습니다.
t1 * d1 + t2 * d2 = p
위의 식을 연립 방정식으로 풀면 t1과 t2를 얻을 수 있습니다.

점 k가 삼각형 내에 있으려면

0 <= t1 <= 1
0 <= t2 <= 1
t1 + t2 <= 1

이 세가지 조건을 만족하면 됩니다.


수학의 세계는 너무나 아름답습니다.  수많은 수학자 여러분께 감사하단 인사를 드립니다. 꾸바닥!