浅谈计算圆周率、黄金比例和自然常数
正文
利用程序计算圆周率、黄金比例和自然常数有很多方法，这里我们举一些小例子管中窥豹。
计算圆周率
1. 统计方法
本质是随机撒点，也叫蒙特·卡罗法（Monte Carlo method）。
源代码：1
2
3
4
5
6
7
8
9
10
11
12
13
14double calc_pi_1(int n)
{
int hits = 0;
for (int i = 0; i < n; i++)
{
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;
if (x * x + y * y < 1.0)
{
hits++;
}
}
return (double)hits / n * 4;
}
运行结果：1
iteration: 1000, result: 3.11200000000000, duration: 32000 ns (Monte Carlo method)
n 取这么小是因为这里只是一个 demo，真正要计算圆周率的话不可能这么小的。
2. 统计方法
如果说蒙特卡罗法本质是随机撒点，那么这个本质就是均匀撒点。
源代码：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15double calc_pi_2(int n)
{
int hits = 0;
for (int x = 0; x * x < n; x++)
{
for (int y = 0; y * y < n; y++)
{
if (x * x + y * y < n)
{
hits++;
}
}
}
return (double)hits / n * 4;
}
运行结果：1
iteration: 1000, result: 3.26000000000000, duration: 1300 ns (Uniform sprinkling method)
虽然结果精确度比蒙特卡罗法稍差，但快了许多，因为蒙特卡罗法每次循环都要调用 rand() 函数来产生伪随机数（不引入特殊设备的话，计算机是无法产生真正的随机数的）。
3. 数学方法
数学是一个强有力的工具，利用级数展开的公式可以快速估算圆周率：
\[ \frac{\pi}{4}=\sum_{n=0}^{\infty}\frac{(-1)^n}{2 n+1} \]
源代码：1
2
3
4
5
6
7
8
9
10double calc_pi_3(int n)
{
double pi = 0, k = 1;
for (int i = 0; i < n; i++)
{
pi += k / (2 * i + 1);
k = -k;
}
return pi * 4;
}
运行结果：1
iteration: 1000, result: 3.14059265383979, duration: 1300 ns (Series summation)
4. 数学方法
上面那个公式效果似乎不尽人意。当然我们可以增加迭代次数，但是否有收敛更快的公式呢？看看下面这个怎么样：
\[ \frac{\pi}{2}=\sum_{n=0}^{\infty}\frac{n!}{(2 n+1)!!} \]
源代码：1
2
3
4
5
6
7
8
9
10double calc_pi_4(int n)
{
double pi = 1, x = 1;
for (int i = 1; i < n; i++)
{
x *= (double)i / (2 * i + 1);
pi += x;
}
return pi * 2;
}
运行结果：1
iteration: 1000, result: 3.14159265358979, duration: 1400 ns (Advanced series summation)
Nice！
可以看到，仅仅一千次迭代，就有如此高的精度，且耗费时间才 1.4 微秒，很不错呢。
计时函数
对了，关于计时，其实很简单：1
2
3
4
5
6
7
8
9
10void time_test(function_t function, const char *message, int n)
{
struct timespec start, end;
clock_gettime(CLOCK_MONOTONIC, &start);
double result = function(n);
clock_gettime(CLOCK_MONOTONIC, &end);
time_t nanoseconds = end.tv_nsec - start.tv_nsec;
printf("iteration: %d, result: %.14lf, duration: %5lld ns (%s)\n", n, result, nanoseconds, message);
}
clock_gettime() 两次调用做差，即得中间过程的运行耗时。
计算黄金比例
1. 数学方法
首先最直接的是利用黄金比例的代数定义来计算：
\[ \frac{a+b}{a} = \frac{a}{b} = \varphi \rightarrow \varphi = \frac{1+\sqrt{5}}{2} = 1.618 033 988 749 ... \]
这个方法并不需要迭代。
源代码：1
2
3
4
5double calc_phi_1(int _)
{
(void)_; // unused
return (1 + sqrt(5)) / 2;
}
运行结果：1
iteration: 1000, result: 1.61803398874989, duration: 100 ns (Algebraic)
2. 数学方法
这里用斐波那契数列（Fibonacci sequence）来计算黄金比例：
\[ t_{n}=1,1,2,3,5,8,13\ldots \]
后一项与前一项的比值：
\[ \frac{1}{1}=1.000, \frac{2}{1}=2.000, \frac{3}{2}=1.500, \frac{5}{3} \approx 1.667, \frac{8}{5}=1.600, \frac{13}{8}=1.625 ... \]
可以看到，随着项数增大，比值越来越接近黄金比例，即：
\[ \lim_{n \rightarrow \infty} \frac{t_{n}}{t_{n-1}}=\varphi \]
源代码：1
2
3
4
5
6
7
8
9
10
11double calc_phi_2(int n)
{
double t1 = 0, t2 = 1;
for (int i = 0; i < n; ++i)
{
double next = t1 + t2;
t1 = t2;
t2 = next;
}
return t2 / t1;
}
运行结果：1
iteration: 1000, result: 1.61803398874989, duration: 1600 ns (Fibonacci sequence)
3. 数学方法
这里要用到连分数（continued fraction）。
第一次看到连分数时是阅读关于拉马努金（Srīṉivāsa Rāmāṉujan）的文章时看见的。此人天赋异禀，堪称鬼才，发现了许多画风奇特的数学公式，比如：
\[ \sqrt{\varphi+2}-\varphi=\frac{e^{-\frac{2 \pi}{5}}}{1+\frac{e^{-2 \pi}}{1+\frac{e^{-4 \pi}}{1+\frac{e^{-6 \pi}}{1+\ddots}}}} \]
把黄金比例 \(\varphi\) 、自然常数 \(e\) 和圆周率 \(\pi\) 联系在了一起。
还有一个计算圆周率收敛极快的公式：
\[ \frac{1}{\pi} = \frac{2 \sqrt{2}}{9801} \sum_{n=0}^{\infty} \frac{(4n)! (1103+26390n)}{(n!)^4 396^{4n}} \]
迭代一次就有 8 位十进制有效数字，迭代两次就达到了 double 精度上限，太牛了……
因为按照公式直接写程序迭代几次之后就会整型上溢，根本没法迭代一千次，所以就不和其他方法进行统一比较了。
还有很多我就不一一展示了，总之那些公式看得我头皮发麻……
连分数就像是一台数学显微镜，它能把放进去的数字的内部构造展现出来，让我们得以一窥精妙的造物之痕。现在，让我们把黄金比例放进去...：
\[ \varphi=1+\frac{1}{1+\frac{1}{1+\frac{1}{1+\ddots}}} \]
化为连分数后，整个式子只有数字 1！
All in one, beautiful!
源代码：1
2
3
4double calc_phi_3(int n)
{
return n == 0 ? 1.618 : 1 + 1 / calc_phi_3(n - 1);
}
一行代码搞定(￣ ▽ ￣)"
运行结果：1
iteration: 1000, result: 1.61803398874989, duration: 25000 ns (Continued fraction)
值得一提的是，递归基例我不知道是多少，一开始设置成 1，结果是对的。我又试了试 0.5，结果还是一样，我有点摸不着头脑了。我再试了一下 0，结果还是一样！有点方，稳住。我把输出的小数位数调到了 14 位，因为 double 有效位数至少 15 位，再把迭代次数调小，这样就能明显地比较每个基例的收敛速度，于是我发现似乎在 1.6 左右时收敛最快……1.618！它的基例，正是它自己本身！这和它的连分数似乎有某种内在的联系，都是“1”……
Σ(っ °Д °;)っ
计算自然常数
1. 数学方法
这里我们用定义直接算：
\[ e=\lim_{n\rightarrow \infty }(1+\frac{1}{n})^n \]
源代码：1
2
3
4double calc_e_1(int n)
{
return pow(1 + 1.0 / n, n);
}
运行结果：1
iteration: 1000, result: 2.71692393223552, duration: 700 ns (Definition)
Oh……似乎精度不够？毕竟这只是定义式，我们接下来换个收敛较快的公式。
2. 数学方法
级数展开，这是一个很优美的公式：
\[ e=\sum_{n=0}^{\infty}{\frac{1}{n!}} \]
源代码：1
2
3
4
5
6
7
8
9
10double calc_e_2(int n)
{
double e = 1, x = 1;
for (int i = 1; i < n; i++)
{
x *= 1.0 / i;
e += x;
}
return e;
}
运行结果：1
iteration: 1000, result: 2.71828182845905, duration: 1700 ns (Series summation)
可以看到，同样的迭代次数，精度明显提高了许多。
3. 统计方法
这种方法是知乎上一位大佬 Milo Yip 在他的文章中写出来的，有点深奥。
源代码：1
2
3
4
5
6
7
8
9
10
11
12
13double calc_e_3(int n)
{
int k = 0;
n /= 2.718; // keep for loop count equals n as far as possible
for (int i = 0; i < n; i++)
{
for (int j = 0; j <= RAND_MAX; j += rand())
{
k++;
}
}
return (double)k / n;
}
运行结果：1
iteration: 1000, result: 2.67029972752044, duration: 16000 ns (Unknown advanced method)
背后的数学原理可以参考：
附录
完整的输出（gcc -O3）：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15Calculate pi:
iteration: 1000, result: 3.11200000000000, duration: 32000 ns (Monte Carlo method)
iteration: 1000, result: 3.26000000000000, duration: 1300 ns (Uniform sprinkling method)
iteration: 1000, result: 3.14059265383979, duration: 1300 ns (Series summation)
iteration: 1000, result: 3.14159265358979, duration: 1400 ns (Advanced series summation)
Calculate phi:
iteration: 1000, result: 1.61803398874989, duration: 100 ns (Algebraic)
iteration: 1000, result: 1.61803398874989, duration: 1600 ns (Fibonacci sequence)
iteration: 1000, result: 1.61803398874989, duration: 25000 ns (Continued fraction)
Calculate e:
iteration: 1000, result: 2.71692393223552, duration: 700 ns (Definition)
iteration: 1000, result: 2.71828182845905, duration: 1700 ns (Series summation)
iteration: 1000, result: 2.67029972752044, duration: 16000 ns (Unknown advanced method)
完整的代码在 GitHub 上： timing_calc_pi_phi_e.c