C/C++ 원주율 계산 프로그램

여기서는 수학 상수 중 하나인 원주율을 계산하는 간단한 프로그램을 소개해볼까 하는데요. 물론 원주율의 값은 이미 잘 알려져 있고 C언어 수학 함수 라이브러리에도 있지만, 단순히 값만 외우는거랑 직접 계산해보는 것은 다르다고 생각했기 때문에 이 포스팅을 하게 되었습니다.

 

반응형

 

원주율을 계산하기 위한 방법으로는 무한 급수나 몬테-카를로 (Monte-Carlo) 방법등 여러가지가 있습니다만, 여기서는 정적분을 통해서 원주율을 구해 보겠습니다.

 

 

이는 역삼각함수 중 하나인 아크탄젠트 함수가 특정한 값을 인자로 받는 경우에, 원주율에 비례하는 값을 돌려준다는 사실을 이용한 방법인데요. 원주율과 삼각함수에 대한 자세한 내용이 궁금하시다면, 다음 포스팅을 참고하면 도움이 되리라 생각합니다.

 

수학 상식 : 원주율과 삼각함수

 

정적분의 값을 계산하기 위해서 사다리꼴 공식 (trapezoidal rule)을 적용하게 될텐데요. 적분구간을 여러개로 나누고, 함수의 형태에 따라 결정되는 여러개의 사다리꼴의 면적에 입각해서 정적분의 근사값을 구하는 방법입니다. 다음 포스팅에 더 자세한 내용이 소개되어 있습니다.

 

사다리꼴 공식을 이용한 C/C++ 수치 적분

 

C언어 소스 코드입니다. 참고로 C++ 에서도 같은 코드를 적용할 수 있습니다.

 

#include<stdio.h>
#include<stdlib.h>
#include<math.h>

double eps_precision = 1e-8;

double func_integrand(double u);

int main(int argc, char *argv[]) {
    double pi_now = 0.;
    double pi_prev;
    int converging = 0;

    unsigned long int nbin_u = 1;
    double u_min = 0.;
    double u_max = 1.;
    double delta_u =
        fabs(u_max - u_min) / (double)nbin_u;

    int istep = 1;
    while (converging == 0) {
        pi_prev = pi_now;

        unsigned int iu;
        if (istep == 1) {
            pi_now = 0.;
            for (iu = 0; iu < nbin_u; iu++) {
                double u0 = u_min + delta_u * (double)iu;
                double u1 = u0 + delta_u;
                pi_now +=
                    0.5 * delta_u * (func_integrand(u0) +
                                     func_integrand(u1));
            }
        } else {
            pi_now = 0.5 * pi_prev;
            for (iu = 0; iu <= nbin_u; iu++) {
                if (iu % 2 == 0) {
                    continue;
                }

                double u_now = u_min + delta_u * (double)iu;
                pi_now +=
                    delta_u * func_integrand(u_now);
            }
        }

        fprintf(stdout,
            "    step %d : pi = %.8f\n", istep, pi_now);
            
        if (fabs(pi_now - pi_prev) <
            0.5 * eps_precision *
            fabs(pi_now + pi_prev)) {
            converging = 1;
        }

        istep += 1;
        nbin_u = 2 * nbin_u;
        delta_u = 0.5 * delta_u;
    }

    fprintf(stdout, "\n");
    fprintf(stdout, "pi from numerical integration\n");
    fprintf(stdout, "  > pi = %.8f\n", pi_now);
    fprintf(stdout, "pi from C math library\n");
    fprintf(stdout, "  > pi = %.8f\n", M_PI);

    return 0;
}

double func_integrand(double u) {
    return 2. / (fabs(u * u) + fabs((1. - u) * (1. - u)));
}

 

실수형 변수 eps_precision을 통해서 원하는 정밀도를 상정하고, 원주율의 값이 그만큼 정확해질때까지 수치적분의 격자 갯수를 증가시키는 방식의 프로그램입니다. 즉 분할된 구간의 갯수를 늘려가면서 리만 합을 계산하고 있죠. 이전 단계와 현재 단계에서의 원주율의 추정치를 저장하기 위한 실수형 변수 pi_prev와 pi_now가 선언되어 있습니다.

 

추가로 수렴을 체크하기 위한 정수형 제어변수 converging이 있는데요. pi_now와 pi_prev의 값을 대조한 다음, 그 차이가 현재 추정치보다 매우 작다면 수렴한다고 가정하고 반복문을 정지시키는 방식입니다. 만약 수렴을 하지 않고 더 정밀한 계산이 필요하다면, 분할된 구간의 갯수 nbin_u를 2배로 늘려서 리만 합을 새로 계산하고 pi_now의 값을 업데이트합니다.

 

코드를 컴파일하고 실행해 보면, 다음과 같은 결과를 얻을 수 있습니다.

 

 

C언어 수학 함수 라이브러리에 정의된 값과 비교했을때, 소수점 7번째 자리까지 일치한다는 것을 볼 수 있습니다.

 

만일 원주율의 값을 더 정밀하게 구하고 싶다면 eps_precision의 값을 줄이면 됩니다만, 그만큼 실행시간은 늘어나겠죠. 뿐만 아니라 리만 합을 위한 분할 구간의 갯수를 저장하는 nbin_u가 unsigned long 정수형 변수이기 때문에, 한계가 존재합니다. 이러한 난점들을 타개하기 위한 방법 중 하나는 병렬 프로그램을 작성하는 것인데요. 정적분의 구간을 여러 개로 나눈 다음 각 프로세스나 스레드에 할당하면, 시간도 절약되고 정밀도도 더 높일 수 있을 것입니다.

 

병렬 프로그래밍을 위한 방법은 대표적으로 OpenMP와 MPI (Massage Passing Interface)가 있습니다. 이들을 활용하여 C언어 또는 C++ 병렬 프로그램을 작성하는 방법은 다음 포스팅들에 소개되어 있습니다.

 

OpenMP 병렬 프로그래밍

OpenMP C/C++ 를 이용한 병렬 프로그래밍

 

MPI 병렬 프로그래밍

MPI C/C++ 를 이용한 병렬 프로그래밍

 

---

 

여기서는 C언어를 통해 프로그램을 작성했지만, C#을 통해서도 원주율을 계산할 수 있습니다. 다음 포스팅에 더 자세한 내용이 소개되어 있습니다.

 

C# 전산수학 - 정적분을 통한 원주율 계산