计算圆周率的C程序

这是很早以前用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