因為對算 pi 的值有點興趣, 所以就寫了一個算看看.
前幾天有人在問開根號至 30 位的, 在這個程式有一段就有用到,
雖然這裡是用整數的, 但是 gmp 中也有浮點數的開平方根,
可以參考看看.
我沒有比較過算出來的值是不是正確的.
另外還有幾個地方可以簡化的, 有興趣的人可以試一試. :)
/*===========================================================================
目的: 用十進位表示 10000 位的 pi, 順便測試 gmp
需求: 先裝 gmp (gnu mp)
方法: 請參考交大資科 BBS 精華區中有關 pi 的算法.
用的是 Chudnovsky 的方法.(discovered by S. Ramanujan)
編譯參數: gcc xxx.c -lgmp
===========================================================================*/
#include <stdio.h>
#include <gmp.h>
int main(void)
{
mpz_t one;
mpz_t pi;
unsigned long k1, k2, k3, k4, k5, k6;
mpz_t t1, t2, t3, t4, t5;
mpz_t up, down;
unsigned long n;
k1=545140134;
k2=13591409;
k3=640320;
k4=100100025;
k5=327843840;
k6=53360;
{
#define DIGITS 10000
char tmp[DIGITS+1];
tmp[0]='1';
for (n=1; n<DIGITS; ++n)
tmp[n]='0';
tmp[n]='\0';
mpz_init_set_str(one, tmp, 10);
}
mpz_init(pi);
mpz_init(t1);
mpz_init(t2);
mpz_init(t3);
mpz_init(t4);
mpz_init(t5);
mpz_init(up); mpz_init(down);
mpz_mul(t1, one, one);
mpz_mul(t2, t1, t1);
mpz_mul_ui(t1, t2, k3);
mpz_sqrt(t2, t1);
mpz_mul_ui(up, t2, k6);
mpz_mul_ui(down, one, k2);
mpz_set(t1, down);
for (n=1; n<DIGITS/14+1; ++n) {
mpz_mul_ui(t2, t1, 6*n-1);
mpz_mul_ui(t3, t2, 6*n-3);
mpz_mul_ui(t2, t3, 6*n-5);
mpz_set_ui(t4, k1);
mpz_mul_ui(t5, t4, n);
mpz_add_ui(t4, t5, k2);
mpz_mul(t3, t2, t4);
mpz_tdiv_q_ui(t1, t3, n);
mpz_tdiv_q_ui(t2, t1, n);
mpz_tdiv_q_ui(t1, t2, n);
mpz_set_ui(t4, k1);
mpz_mul_ui(t5, t4, n-1);
mpz_add_ui(t4, t5, k2);
mpz_tdiv_q(t2, t1, t4);
mpz_tdiv_q_ui(t3, t2, k4);
mpz_tdiv_q_ui(t2, t3, k5);
mpz_neg(t1, t2);
mpz_add(t2, t1, down);
mpz_set(down, t2);
}
mpz_tdiv_q(pi, up, down);
mpz_out_str(stdout, 10, pi);
printf("\n");
mpz_clear(t1);
mpz_clear(t2);
mpz_clear(t3);
mpz_clear(t4);
mpz_clear(t5);
mpz_clear(pi);
mpz_clear(up);
mpz_clear(down);
mpz_clear(one);
}
--
※ Origin: 程式設計樂園 ◆ From: g863304.CHING.ab.nthu.edu.tw