一、公式

二、代码

[cpp]
#include <iostream>
#include <stdlib.h>
#define TAIL 10

using namespace std;

void output(int *num, int len) {
cout << num[0] << ‘.’ << endl;
for (int i = 1; i < len; i++) {
cout << num[i];
if (!(i % 50)) cout << " : " << i << endl;
if (!(i % 500)) cout << endl;
if (!(i % 10) && (i % 50)) cout << ‘ ‘;
}
cout << endl;
}

int main(int argc, char *argv[]) {
int a = 1, len, reg, carry;
int decimal;
int *t, *pi;
int eps = 0; // precision

if (argc != 2) {
cout << "Usage: " << argv[0] << " num_of_digits" << endl;
exit(EXIT_FAILURE);
}
decimal = atoi(argv[1]);
len = decimal + TAIL;

t = new int[len];
pi = new int[len];

cout << "Now we start calculate the value of Pi, please wait…" << endl;

for (int i = 0; i < len; i++) {
pi[i] = 0;
t[i] = 0;
}

pi[0] = 2;
t[0] = 2;

while (eps < len – 1) {
carry = 0;
for (int i = len – 1; i >= 0; i–) // t = t * a;
{
reg = t[i] * a + carry;
t[i] = reg % 10;
carry = reg / 10;
}

carry = 0;
for (int i = 0; i < len; i++) // t = t / (2a + 1)
{
reg = t[i] + carry * 10;
t[i] = reg / (2 * a + 1);
carry = reg % (2 * a + 1);
}

carry = 0;
for (int i = len – 1; i >= 0; i–) // pi = pi + t;
{
reg = pi[i] + t[i] + carry;
pi[i] = reg % 10;
carry = reg / 10;
}

// check precision
for (int i = 1; i < len; i++) {
if (t[i]) {
eps = i;
break;
}
}
a++;

}

output(pi, decimal);

delete[] t;
delete[] pi;
return 0;
}

[/cpp]

三、测试

$ time ./pi 1000
Now we start calculate the value of Pi, please wait...
Pi = 3.
1415926535 8979323846 2643383279 5028841971 6939937510 : 50
5820974944 5923078164 0628620899 8628034825 3421170679 : 100
8214808651 3282306647 0938446095 5058223172 5359408128 : 150
4811174502 8410270193 8521105559 6446229489 5493038196 : 200
4428810975 6659334461 2847564823 3786783165 2712019091 : 250
4564856692 3460348610 4543266482 1339360726 0249141273 : 300
7245870066 0631558817 4881520920 9628292540 9171536436 : 350
7892590360 0113305305 4882046652 1384146951 9415116094 : 400
3305727036 5759591953 0921861173 8193261179 3105118548 : 450
0744623799 6274956735 1885752724 8912279381 8301194912 : 500

9833673362 4406566430 8602139494 6395224737 1907021798 : 550
6094370277 0539217176 2931767523 8467481846 7669405132 : 600
0005681271 4526356082 7785771342 7577896091 7363717872 : 650
1468440901 2249534301 4654958537 1050792279 6892589235 : 700
4201995611 2129021960 8640344181 5981362977 4771309960 : 750
5187072113 4999999837 2978049951 0597317328 1609631859 : 800
5024459455 3469083026 4252230825 3344685035 2619311881 : 850
7101000313 7838752886 5875332083 8142061717 7669147303 : 900
5982534904 2875546873 1159562863 8823537875 9375195778 : 950
1857780532 1712268066 1300192787 6611195909 216420198

real	0m5.230s
user	0m0.184s  //用户空间的CPU时间
sys	0m0.000s

四、分析

1.逼近过程

下表展示了临时变量t以及待求变量pi随自变量a的迭代变化过程。
可见，随着t越来越小，pi也越来越精确。

2.时间复杂度

由以下图表可以看出，时间复杂度为O(n^2)
n<1000

n>1000