/*
名稱: pi.c
作者: 洪朝貴 http://www.cyut.edu.tw/~ckhung/
功能: 從命令列上讀入 prec, 計算 pi=3.1416 ... 精準至 prec 個 bits.
需求: 使用到 GNU Mulitiple Precision (gmp) 程式庫,
編譯: gcc -Wall pi.c -lgmp
*/

#include <stdio.h>
#include <gmp.h>				/* 要在 include stdio.h 之後 */
#include <assert.h>
#include <stdlib.h>

int main(int argc, char * argv[])
/* 演算法: 利用 pi ~= 2^k * sin (pi / 2^k) 所導出的遞迴關係式: 來計算 pi */
{				/* pi(k+1) = pi(k) / cos(pi / 2^(k+1))   */
    int step, prec;             /* cos(th/2) = sqrt((1 + cos(th) / 2)    */
    mpf_t c, pi, old_pi;

    assert(2 >= argc);
    prec = 64;
    if (2 == argc) {
	prec = atoi(argv[1]);
	assert(prec > 0);
    }
    mpf_set_default_prec(prec);			/* 設定內定的精準度為 prec */

    mpf_init(pi);				/* 起始高精準度變數 */
    mpf_init(old_pi);
    mpf_init(c);
    mpf_set_ui(pi, 2);				/* pi 的近似值 */
    mpf_set_ui(old_pi, 0);			/* pi 上次的近似值 */
    mpf_set_ui(c, 0);				/* cosine 值 */
    for (step=0; mpf_cmp(old_pi,pi); ++step) {	/* 只要 pi 的值還在改變... */
	mpf_set(old_pi, pi);			/* old_pi <-- pi */
	mpf_add_ui(c, c, 1);
	mpf_div_ui(c, c, 2);
	mpf_sqrt(c, c);				/* c <-- sqrt((c+1)/2) */
	mpf_div(pi, pi, c);			/* pi <-- pi / c; */
	printf("%3d: ", step);
	mpf_out_str(stdout, 10, 0, pi);		/* 印出 pi */
	putchar('\n');
    }
    mpf_clear(c);				/* 釋放高精準度變數的空間 */
    mpf_clear(old_pi);
    mpf_clear(pi);
    return 0;
}