这是很早以前用C语言写的一个计算圆周率的程序, 算法是用泰勒公式计算反正切值。在命令行不跟参数执行该程序则使用Gauss公式计算前1000位圆周率的值,如果带一个命令行参数,则该值为要计算的位数。如果还有第二个命令行参数,则使用Stomer公式计算,可作为验算。因为该程序只涉及到纯数学计算,可以在Linux、Unix、Windows等操作系统下编译并运行。当时写这个程序时,int是2个字节的,现在大多数的C编译器int都是4个字节,不过这不影响程序的正确性。
#include
#include
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])) 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 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 |