這是很早以前用C語言寫的一個計算圓周率的程序, 算法是用泰勒公式計算反正切值。在命令行不跟參數執行該程序則使用Gauss公式計算前1000位圓周率的值，如果帶一個命令行參數，則該值為要計算的位數。如果還有第二個命令行參數，則使用Stomer公式計算，可作為驗算。因為該程序只涉及到純數學計算，可以在Linux、Unix、Windows等操作系統下編譯並運行。當時寫這個程序時，int是2個字節的，現在大多數的C編譯器int都是4個字節，不過這不影響程序的正確性。
#include <stdio.h>
#include <stdlib.h>
main(int argc, char * argv[])
{
long * pi, * t, m, n, r, s;
int t0[][3] = {48, 32, 20, 24, 8, 4}, k0[][3] = {1, 1, 0, 1, 1, 1};
int n0[][3] = {18, 57, 239, 8, 57, 239}, d, i, j, k, p, q;
d = (argc > 1) ? (((i = atoi(argv[1])) < 0) ? 0 : i) : 1000;
q = (argc > 2) ? 1 : 0;
printf("%s\n\n", "Nature (R) Pi value compute Program (C) Tue 1999.11.30");
printf("pi= %s%d * arctg(1/%d) %s %d * arctg(1/%d) %s %d * arctg(1/%d) [%s]\n",
k0[q][0] ? "" : "-", t0[q][0], n0[q][0], k0[q][1] ? "+" : "-", t0[q][1],
n0[q][1], k0[q][2] ? "+" : "-", t0[q][2], n0[q][2], q ? "Stomer" : "Gauss");
if ((t = (long *)calloc((d += 5) + 1, sizeof(long))) == NULL) return 1;
if ((pi = (long *)calloc(d + 1, sizeof(long))) == NULL) return 2;
for (i = d; i >= 0; i--) pi[i] = 0;
for (p = 0; p < 3; p++) {
for (k=k0[q][p], n=n0[q][p], t[i=j=d]=t0[q][p], i--; i >= 0; i--) t[i] = 0;
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % n;
t[i] = m / n;
k ? (pi[i] += t[i]) : (pi[i] -= t[i]);
}
while (j > 0 && t[j] == 0) j--;
for (k = !k, s = 3, n *= n; j > 0; k = !k, s += 2) {
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % n;
t[i] = m / n;
}
while (j > 0 && t[j] == 0) j--;
for (r = 0, i = j; i >= 0; i--) {
r = (m = 10 * r + t[i]) % s;
m /= s;
k ? (pi[i] += m) : (pi[i] -= m);
}
}
}
for (n = i = 0; i <= d; pi[i++] = r) {
n = (m = pi[i] + n) / 10;
if ((r = m % 10) < 0) r += 10, n--;
}
printf("pi= %ld.", pi[d]);
for (i = d - 1; i >= 5; i--)
printf("%ld%s", pi[i], ((m = d - i + 5) % 65) ? ((m % 5) ? "" : " ") : "\n");
printf("%sDIGITS: %d\n", (m % 65) ? "\n" : "", d - 5);
return 0;
}
下面是程序運行結果：
Nature (R) Pi value compute Program (C) Tue 1999.11.30
pi= 48 * arctg(1/18) + 32 * arctg(1/57) - 20 * arctg(1/239) [Gauss]
pi= 3.14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 58209 74944
59230 78164 06286 20899 86280 34825 34211 70679 82148 08651 32823 06647 09384
46095 50582 23172 53594 08128 48111 74502 84102 70193 85211 05559 64462 29489
54930 38196 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 45648
56692 34603 48610 45432 66482 13393 60726 02491 41273 72458 70066 06315 58817
48815 20920 96282 92540 91715 36436 78925 90360 01133 05305 48820 46652 13841
46951 94151 16094 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548
07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 98336 73362 44065
66430 86021 39494 63952 24737 19070 21798 60943 70277 05392 17176 29317 67523
84674 81846 76694 05132 00056 81271 45263 56082 77857 71342 75778 96091 73637
17872 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 42019 95611
21290 21960 86403 44181 59813 62977 47713 09960 51870 72113 49999 99837 29780
49951 05973 17328 16096 31859 50244 59455 34690 83026 42522 30825 33446 85035
26193 11881 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 59825
34904 28755 46873 11595 62863 88235 37875 93751 95778 18577 80532 17122 68066
13001 92787 66111 95909 21642 01989
DIGITS: 1000
