#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <gmp.h>
#include <time.h>

#define MPF_PREC    16384       /* 仮数部精度、単位は bit */
#define MPF_DIGITS  4931        /* mpf_get_str 桁数 */
#define MPF_LIMIT   "1e-4931"   /* atan 計算ループの終了条件 */

/*
    m_atan = mi * atan(1/ni)　その１
    i が偶数のとき atan()、i が奇数のとき -atan()
*/
void em_atan(mpf_t m_atan, unsigned int mi, unsigned int ni, unsigned int i)  
{
    clock_t t0, t1;                /* 時間計測用 */
    unsigned int n2i = ni * ni;    /* n2 の計算の簡略化に使用します */
    mpf_t n, n2, ta, limit;
    mpf_init_set_ui(n, 1);         /* n = 1, 3, 5, ... */
    mpf_init_set_ui(n2, ni);       /* n2 = ni, ni**3, ni**5, ... */
    mpf_init(ta);
    mpf_init_set_str(limit, MPF_LIMIT, 10);

    t0 = clock();
    for (; ; i += 1) {                  /* i = [ 0, ] 1, 2, ... */
        mpf_mul(ta, n, n2);             /* ta = n * n2 */
        mpf_ui_div(ta, mi, ta);         /* ta = mi / ta */
        if (i & 1)                      /* i が奇数のときは減算 */
            mpf_sub(m_atan, m_atan, ta);
        else
            mpf_add(m_atan, m_atan, ta);
        if (mpf_cmp(ta, limit) < 0)     /* 終了判断 */
            break;
        mpf_add_ui(n, n, 2);            /* n += 2 */
        mpf_mul_ui(n2, n2, n2i);        /* n2 *= n2i */
    }
    t1 = clock();
    printf("%gms\n", (double)1000.0 * (t1 - t0) / CLOCKS_PER_SEC); 
    mpf_clear(n);
    mpf_clear(n2);
    mpf_clear(ta);
    mpf_clear(limit);
}

/*
    m_atan = mi * atan(1/ni)　その２
    i が正のとき atan()、そうでないとき -atan()
*/
void ef_atan(mpf_t m_atan, unsigned int mi, unsigned int ni, int i)  
{
    clock_t t0, t1;               /* 時間計測用 */
    unsigned int d1 = 1;          /* d1 = 1, 3, 5, 7, ... */
    unsigned int n2i = ni * ni;   /* n2i = ni ** 2 */

    mpf_t n1, limit, tn1, tn2;
    mpf_init_set_ui(n1, ni);      /* n1 = ni**1, ni**3, ni**5, ... */
    mpf_init_set_str(limit, MPF_LIMIT, 10);
    mpf_init(tn1);
    mpf_init(tn2);

    t0 = clock();
    if (i > 0) {
        for (;;) {
            mpf_mul_ui(tn1, n1, d1);
            mpf_ui_div(tn2, mi, tn1);      /* tn2 = mi / (n1 * d1) */
            mpf_add(m_atan, m_atan, tn2);  /* m_atan += tn2 */

            d1 += 2;
            mpf_mul_ui(tn1, n1, n2i);      /* 乗除算は in-place 演算を避ける */
            mpf_set(n1, tn1);              /* n1 *= n2i */

            mpf_mul_ui(tn1, n1, d1);
            mpf_ui_div(tn2, mi, tn1);      /* tn2 = mi / (n1 * d1) */
            mpf_sub(m_atan, m_atan, tn2);  /* m_atan -= tn2 */

            if (mpf_cmp(tn2, limit) < 0)   /* 終了判断 */
                break;

            d1 += 2;
            mpf_mul_ui(tn1, n1, n2i);
            mpf_set(n1, tn1);              /* n1 *= n2i */
        }
    } else {
        for (;;) {
            mpf_mul_ui(tn1, n1, d1);
            mpf_ui_div(tn2, mi, tn1);      /* tn2 = mi / (n1 * d1) */
            mpf_sub(m_atan, m_atan, tn2);  /* m_atan -= tn2 */

            d1 += 2;
            mpf_mul_ui(tn1, n1, n2i);      /* 乗除算は in-place 演算を避ける */
            mpf_set(n1, tn1);              /* n1 *= n2i */

            mpf_mul_ui(tn1, n1, d1);
            mpf_ui_div(tn2, mi, tn1);      /* tn2 = mi / (n1 * d1) */
            mpf_add(m_atan, m_atan, tn2);  /* m_atan += tn2 */

            if (mpf_cmp(tn2, limit) < 0)   /* 終了判断 */
                break;

            d1 += 2;
            mpf_mul_ui(tn1, n1, n2i);
            mpf_set(n1, tn1);              /* n1 *= n2i */
        }
    }
    t1 = clock();
    printf("%gms\n", (double)1000.0 * (t1 - t0) / CLOCKS_PER_SEC);
    mpf_clear(n1);
    mpf_clear(limit);
    mpf_clear(tn1);
    mpf_clear(tn2);
}

/*
    数値文字列を見やすく表示
*/
void disp_mstr(char* s)
{
    int i, j, k;
    char cb[12], *d;
    s += 1;                    /* 先頭の '3' をスキップする */
    for (i = 0; *s; i += 1) {
         d = cb;
         for (j = 0; j < 10 && *s; j += 1)
             *d++ = *s++;
         strcpy(d, " ");
         k = i / 10;
         if (k >= 0) {         /* 比較する値を大きくすると、後ろのほうのみ表示できます */
             printf(cb);
             if (i % 10 == 9)
                 printf(": %d\n", k + 1);
         }
    }
    printf("\n");
}

/*
    マチンの公式を使ってπの近似値を求めます
*/
int pi_mpf()
{
    void* sa = malloc(5000);    /* 出力文字列のスペース */
    clock_t t0, t1;             /* 時間計測用 */
    mp_exp_t ve;

    mpf_t m_atan;
    mpf_set_default_prec(MPF_PREC); /* 単位は bit */
    mpf_init(m_atan);               /* m_atan = 0 */

/* 方法その１ */
    em_atan(m_atan, 16, 5, 0);      /* +16 * atan(1/5) */
    em_atan(m_atan, 4, 239, 1);     /* -4 * atan(1/239) */

    t0 = clock();
    char* s = mpf_get_str(sa, &ve, 10, 48, m_atan);  /* 48 桁と指定 */
    t1 = clock();
    printf("%gms\n", (double)1000.0 * (t1 - t0) / CLOCKS_PER_SEC); 
    if (s)
        disp_mstr(s);

/* 方法その２ */
    mpf_set_d(m_atan, 0);        /* m_atan = 0 */
    ef_atan(m_atan, 16, 5, 1);
    ef_atan(m_atan, 4, 239, -1);
    s = mpf_get_str(sa, &ve, 10, MPF_DIGITS, m_atan);
    if (s)
        disp_mstr(s);

    mpf_clear(m_atan);
    if (sa)
        free(sa);
    return 0;
}

int main(int argc, char* argv[])
{
    pi_mpf();
    return 0;
}