#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
// MPI 헤더 파일
#include<mpi.h>

int n_size_;  // 프로세스의 총 갯수
int n_rank_;  // 각 프로세스에 부여된 랭크
// 정밀도
double eps_precision = 1e-12;

// 적분 대상이 되는 함수
double func_integrand(double u);

int main(int argc, char *argv[]) {
    // MPI 초기화
    MPI_Init(&argc, &argv);

    // 프로세스 총 갯수 및 각 프로세스의 랭크
    MPI_Comm_size(MPI_COMM_WORLD, &n_size_);
    MPI_Comm_rank(MPI_COMM_WORLD, &n_rank_);

    if (n_rank_ == 0) {
        // 랭크가 0 인 프로세스
        fprintf(stdout,
            "We have %d processes.\n", n_size_);
        fprintf(stdout, "\n");
    }

    // 모든 프로세스가 여기에 도달할 때 까지 대기
    MPI_Barrier(MPI_COMM_WORLD);

    // MPI 통신을 위한 변수들
    int tag;
    MPI_Status status;

    // 현재 단계의 원주율 값
    double pi_now = 0.;
    // 이전 단계의 원주율 값
    double pi_prev;
    // 현재 단계의 각 프로세스의 기여분
    double pi_rank = 0.;
    // 이전 단계의 각 프로세스의 기여분
    double pi_rank_prev;
    // 수렴 여부를 체크하기 위한 제어 변수
    int converging = 0;

    /* 각 프로세스에서 수치적분을 위한 구간 및
     * 격자 간격 */
    unsigned long int nbin_u = 1;
    double u_min = (double)n_rank_ / (double)n_size_;
    double u_max = u_min + 1. / (double)n_size_;
    double delta_u = fabs(u_max - u_min);

    int istep = 1;
    /* 지정한 정밀도 이내에서 수렴할 때 까지
     * 반복문 실행 */
    while (converging == 0) {
        pi_prev = pi_now;

        pi_rank_prev = pi_rank;
        pi_rank = 0.;
        unsigned int iu;
        // 수치적분 계산
        if (istep == 1) {
            for (iu = 0; iu < nbin_u; iu++) {
                double u0 = u_min + delta_u * (double)iu;
                double u1 = u0 + delta_u;
                pi_rank +=
                    0.5 * delta_u * (func_integrand(u0) +
                                     func_integrand(u1));
            }
        } else {
            pi_rank = 0.5 * pi_rank_prev;
            for (iu = 0; iu <= nbin_u; iu++) {
                if (iu % 2 == 0) {
                    continue;
                }

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

        pi_now = 0.;
        if (n_rank_ == 0) {
            // 랭크가 0 인 프로세스

            /* 모든 프로세스의 기여분들을 취합하여
             * 원주율의 값 계산 */
            pi_now = pi_rank;
            for (int irank = 1; irank < n_size_; irank++) {
                tag = 1000 + irank;
                double pi_add;
                MPI_Recv(&pi_add, 1, MPI_DOUBLE, irank,
                         tag, MPI_COMM_WORLD, &status);
                pi_now += pi_add;
            }
            fprintf(stdout,
                "    step %d : pi = %.12f\n", istep, pi_now);
        } else {
            // 랭크가 0 이 아닌 프로세스
            tag = 1000 + n_rank_;
            MPI_Send(&pi_rank, 1, MPI_DOUBLE, 0,
                     tag, MPI_COMM_WORLD);
        }

        if (n_rank_ == 0) {
            // 랭크가 0 인 프로세스

            // 수렴 체크
            if (fabs(pi_now - pi_prev) <
                    0.5 * eps_precision *
                    fabs(pi_now + pi_prev)) {
                converging = 1;
            }

            for (int irank = 1; irank < n_size_; irank++) {
                tag = 2000 + irank;
                MPI_Send(&converging, 1, MPI_INT, irank,
                         tag, MPI_COMM_WORLD);
            }
        } else {
            // 랭크가 0 이 아닌 프로세스
            tag = 2000 + n_rank_;
            MPI_Recv(&converging, 1, MPI_INT, 0,
                     tag, MPI_COMM_WORLD, &status);
        }

        istep += 1;
        // 적분 구간의 격자 갯수를 2배로 증가
        nbin_u = 2 * nbin_u;
        delta_u = 0.5 * delta_u;
    }

    // 모든 프로세스가 여기에 도달할 때 까지 대기
    MPI_Barrier(MPI_COMM_WORLD);

    if (n_rank_ == 0) {
        // 랭크가 0 인 프로세스
        fprintf(stdout, "\n");
        fprintf(stdout, "pi from numerical integration\n");
        fprintf(stdout, "  > pi = %.12f\n", pi_now);
        fprintf(stdout, "pi from C math library\n");
        fprintf(stdout, "  > pi = %.12f\n", M_PI);
    }

    // MPI 종료
    MPI_Finalize();

    return 0;
}

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