/*
    mpf_t と double 演算速度の比較
    仮数部精度 128 bit の場合
*/

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

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

/*
    m_atan = mi * atan(1/ni)
    i が正のとき atan()、そうでないとき -atan()
*/
void mf_atan(mpf_t m_atan, unsigned int mi, unsigned int ni, int i)  
{
    int j;                        /* 50000 回ループ用 */
    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, sav_matan;
    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);
    mpf_init(sav_matan);
    mpf_set(sav_matan, m_atan);   /* 値をセーブしておく */

    t0 = clock();
    for (j = 0; j < 50000; j += 1) {
        d1 = 1;                       /* 変数の値を復帰する */
        mpf_set(m_atan, sav_matan);
        mpf_set_ui(n1, ni);
        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(n1, n1, n2i);       /* 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(n1, n1, n2i);       /* 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(n1, n1, n2i);       /* 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(n1, n1, n2i);       /* n1 *= n2i */
            }
        }
    }
    t1 = clock();
    printf("%3u: %gms\n", ni, (double)1000.0 * (t1 - t0) / CLOCKS_PER_SEC);
    mpf_clear(n1);
    mpf_clear(limit);
    mpf_clear(tn1);
    mpf_clear(tn2);
    mpf_clear(sav_matan);
}

/*
    mf_atan() の double 版
*/
void df_atan(double* d_atan, unsigned int mi, unsigned int ni, int i)  
{
    int j;                        /* 17 * 50000 回ループ用 */
    clock_t t0, t1;               /* 時間計測用 */
    unsigned int d1 = 1;          /* d1 = 1, 3, 5, 7, ... */
    unsigned int n2i = ni * ni;   /* n2i = ni ** 2 */

    double n1, limit, tn2, loc_atan;
    n1 = ni;                      /* n1 = ni**1, ni**3, ni**5, ... */
    limit = 1e-16;
    t0 = clock();
    for (j = 0; j < 17 * 50000; j += 1) {
        loc_atan = *d_atan;       /* 変数の値を復帰する */
        d1 = 1;
        n1 = ni;
        if (i > 0) {
            for (;;) {
                loc_atan += mi / (n1 * d1);
                d1 += 2;
                n1 *= n2i;

                tn2 = mi / (n1 * d1);
                loc_atan -= tn2;
                if (tn2 < limit)   /* 終了判断 */
                    break;
                d1 += 2;
                n1 *= n2i;
            }
        } else {
            for (;;) {
                loc_atan -= mi / (n1 * d1);
                d1 += 2;
                n1 *= n2i;

                tn2 = mi / (n1 * d1);
                loc_atan += tn2;
                if (tn2 < limit)   /* 終了判断 */
                    break;
                d1 += 2;
                n1 *= n2i;
            }
        }
    }
    *d_atan = loc_atan;     /* 計算結果 */
    t1 = clock();
    printf("%3u: %gms\n", ni, (double)1000.0 * (t1 - t0) / CLOCKS_PER_SEC);
}

/*
    数値文字列を見やすく表示
*/
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");
}

/*
    πを求めるマチンの atan() 展開式を使って、mpf_t と
    double 型での演算速度を比較します
*/
void pi_mpf_d()
{
    mp_exp_t ve;
    mpf_t m_atan;
    void* sa = malloc(5000);        /* mpf_t 出力文字列のスペース */
    mpf_set_default_prec(MPF_PREC); /* 単位は bit */
    mpf_init(m_atan);               /* m_atan = 0 */

    printf("mpf_t\n");
    mf_atan(m_atan, 16, 5, 1);
    mf_atan(m_atan, 4, 239, -1);
    char* s = mpf_get_str(sa, &ve, 10, MPF_DIGITS, m_atan);
    if (s)
        disp_mstr(s);
    mpf_clear(m_atan);
    if (sa)
        free(sa);

    printf("\ndouble\n");
    double d_atan = 0;               /* double 型での演算 */
    df_atan(&d_atan, 16, 5, 1);
    df_atan(&d_atan, 4, 239, -1);
    printf("%.15e\n", d_atan);
}

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