前言

数学中的圆周率 π 我们在很小的时候，就学习了。

圆周率 = 圆的周长 / 圆的直径。

我们也都背的很熟，π = 3.1415926…

那么今天我们来进一步学习一下，如果想尽可能精确地计算圆周率 pi，应该怎么实现呢？

历史与发展

回头看中国

圆周率（π）的历史在中国可以追溯到古代，当时的数学家们已经开始探索并计算这个重要的数学常数。

以下是中国圆周率历史的一些重要里程碑：

早期记载与探索

《周髀算经》：这是一部最晚成书于公元前1世纪的中国早期数学著作，其中明确指出：“圆径一而周三”，即π的近似值为3。虽然这个比值不够精确，但它包含了明确的圆周率概念。
《九章算术》：成书于东汉时期，是中国第一部数学专著。在计算圆的面积、圆柱和圆锥的体积时，使用了π=3的近似值。

祖冲之的贡献

祖冲之（约公元429年—公元500年）：南北朝时期的著名数学家和天文学家，他在数学上的重要贡献之一是计算出圆周率π的值在3.1415926和3.1415927之间。祖冲之提出了π的约率22/7和密率355/113，其中密率值比欧洲早了1000多年。
祖冲之的计算方法：虽然具体的方法没有详细记载，但祖冲之的成就表明他使用了非常先进的数学技巧。他的密率355/113，是一个精确到小数点后七位的近似值，这一成就在当时是无与伦比的。

后世的影响与纪念

月球上的祖冲之山：为了纪念祖冲之的贡献，月球上的一座环形山被命名为“祖冲之山”。
学术研究与讨论：中国科学家茅以升在20世纪初对圆周率的历史进行了研究，并在1917年发表文章《中国圆周率略史》，提出了将圆周率称为“祖率”的建议。
国际认可：日本数学家三上义夫在《中日数学发展史》中建议将355/113称为“π的祖冲之分数值”，以表彰祖冲之的贡献。

抬头看世界

圆周率（π）是一个数学常数，表示圆的周长与直径的比例。

它是一个无理数，意味着它不能表示为两个整数的比例，其小数部分是无限不循环的。

1. 阿基米德方法（几何方法）

阿基米德是第一个使用几何方法来计算圆周率近似值的数学家。

他采用了一种称为“穷竭法”的技巧，通过计算多边形的周长来逼近圆的周长。阿基米德使用正多边形（例如正六边形、正十二边形等）来逼近圆，边数越多，逼近的精度越高。

2. 蒙特卡洛方法（统计方法）

蒙特卡洛方法是一种统计模拟方法，可以用来近似计算圆周率。该方法涉及随机抛洒点到一个正方形内，并计算落在内切圆内的点的比例。通过足够多次的试验，可以得到圆周率的近似值。这个方法的基本原理是基于几何概率。

3. 莱布尼茨级数（解析方法）

莱布尼茨在17世纪提出了一个著名的无穷级数，可以用来计算圆周率：

π = 4 * (1 - 1/3 + 1/5 - 1/7 + 1/9 - …)

这个级数是交错级数，其中每一项的分母是奇数，符号交替为正负。这个级数收敛较慢，需要计算很多项才能得到π的较高精度的近似值。

4. 楚德诺夫斯基算法（快速收敛算法）

楚德诺夫斯基算法是20世纪发展的一种快速收敛的算法，用于计算圆周率。这个算法基于一个特定的级数，可以快速收敛到π的准确值。它是目前已知最快的算法之一，用于计算圆周率的高精度数值。

这些方法只是计算圆周率的众多方法中的一部分。随着计算技术的发展，现在已经有更多高效和精确的方法来计算圆周率，甚至可以计算到数十亿甚至更多位数的小数。

阿基米德方法（几何方法）

阿基米德是古希腊的一位伟大数学家、物理学家和工程师，他对圆周率π的研究做出了重要贡献。

阿基米德方法，也称为穷竭法（Exhaustion Method），是一种几何方法，用于计算圆的面积和体积以及其他几何形状的体积。

这种方法通过无限逼近的方式，精确计算出圆周率的值。

穷解法

阿基米德的穷竭法

阿基米德在其著作《论球与圆柱》中详细描述了穷竭法。他的目标是计算圆的面积和球体的体积。

为了做到这一点，他设想了一种方法，通过构造正多边形来逼近圆形和球形。

正多边形的构造：阿基米德首先构造了一个正六边形，使其与一个给定的圆相切（内切）。这意味着正六边形的所有顶点都位于圆上。
逼近圆的周长：通过增加正多边形的边数，阿基米德能够更精确地逼近圆的周长。随着边数的增加，正多边形的边越来越短，形状越来越接近圆形。理论上，当边数无限增加时，正多边形将与圆完全重合。
计算圆的面积：阿基米德通过计算正多边形的面积来逼近圆的面积。正多边形的面积可以通过其边长和中心角来计算。通过不断增加边数，正多边形的面积也越来越接近圆的面积。
圆周率的计算：阿基米德使用这种方法，计算出了圆周率的一个界限。他发现，圆的周长与直径的比值（即圆周率π）必须大于某个值，同时小于另一个值。通过这种方法，他得到了π的一个近似值，这个值在3.14和3.142之间。

蒙特卡洛方法（统计方法）

说明

蒙特卡洛方法是一种基于随机抽样来解决数学和物理问题的方法，特别是在计算复杂或无法解析求解的问题时非常有效。

这种方法得名于摩纳哥的蒙特卡洛赌场，因为它的核心思想与赌博中的随机性有关。

在计算圆周率π的值时，蒙特卡洛方法提供了一种独特的途径。

蒙特卡洛方法的基本原理

蒙特卡洛方法依赖于随机数和概率统计原理。它通过生成大量随机样本来估计一个量的数值。在计算π的情况下，这个方法通常涉及以下步骤：

选择一个单位正方形：首先，我们在平面上画一个边长为1的正方形。这个正方形的面积是1。
内切圆：在正方形内部画一个内切圆，即圆与正方形的每条边都相切。对于单位正方形，内切圆的半径是1/√2，因此圆的面积是π/4。
随机点生成：接下来，我们在单位正方形内生成大量随机点。这些点的坐标（x, y）是在0到1之间随机选择的。
点与圆的关系判断：对于每个随机点，我们检查它是否落在内切圆内。如果点的坐标满足x^2 + y^2 ≤ 1的条件，那么这个点就在圆内。
概率和面积估计：由于正方形和内切圆的面积已知，我们可以计算圆内点的比例。这个比例等于圆的面积除以正方形的面积。因此，我们可以用这个比例来估计π的值。

蒙特卡洛方法的实现

在实践中，蒙特卡洛方法的实现通常涉及以下步骤：

使用随机数生成器产生一系列随机点的坐标。
对于每个点，检查它是否在内切圆内。
计算落在圆内的点数。
使用圆内点数与总点数的比值乘以4，得到π的近似值。

蒙特卡洛方法的优点

蒙特卡洛方法的优点在于它的简单性和灵活性。它不需要复杂的数学公式或解析解，只需要能够生成随机数的计算机程序。此外，蒙特卡洛方法可以应用于各种问题，不仅限于计算π，还包括求解积分、优化问题、模拟物理过程等。

蒙特卡洛方法的局限性

尽管蒙特卡洛方法非常强大，但它也有一些局限性。最主要的是，蒙特卡洛方法是一种统计方法，它的结果是基于概率的，因此存在一定的误差范围。此外，蒙特卡洛方法的精度依赖于样本数量，需要足够多的随机点才能得到高精度的结果。

总的来说，蒙特卡洛方法是一种非常有用的工具，它通过随机抽样和概率统计原理，为我们提供了一种解决复杂问题的有效途径。

在计算圆周率π的值时，蒙特卡洛方法展示了其独特的优势和实用性。

程序代码

我们通过简单的程序实现模拟一下：

package com.github.houbb.how.calc.pi;

import java.util.Random;

// 作者：老马啸西风
public class MonteCarloPiEstimation {

    public static void main(String[] args) {
        int totalPoints = 10000000; // 总随机点数
        Random random = new Random();
        int pointsInsideCircle = 0; // 圆内的点数

        // 生成随机点并统计落在圆内的点数
        for (int i = 0; i < totalPoints; i++) {
            double x = random.nextDouble(); // 生成0到1之间的随机x坐标
            double y = random.nextDouble(); // 生成0到1之间的随机y坐标

            // 检查点是否在单位圆内 (x^2 + y^2 <= 1)
            if (x * x + y * y <= 1) {
                pointsInsideCircle++;
            }
        }

        // 计算π的近似值
        double piEstimate = 4.0 * pointsInsideCircle / totalPoints;
        System.out.println("Estimated value of Pi: " + piEstimate);
    }

}

1000W 次模拟，近似值为

Estimated value of Pi: 3.1414572

莱布尼茨级数

说明

莱布尼茨级数是一个著名的交错级数，它以德国数学家和哲学家戈特弗里德·威廉·莱布尼茨的名字命名。这个级数是π的一个无穷级数表示，具有重要的数学和历史意义。

莱布尼茨级数的表达式如下：

pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - …

这个级数是一个交错级数，其中每一项的分母是奇数，符号交替为正负。莱布尼茨级数可以用来计算π的近似值，但由于它收敛得相对较慢，因此在实际计算中通常只用于教学和演示目的。

莱布尼茨级数的收敛性

莱布尼茨级数的收敛性可以通过比较测试来证明。尽管收敛速度较慢，但它确实收敛到π/4。每一项的大小随着分母的增大而减小，且符号交替变化。这意味着，随着级数项数的增加，总和逐渐接近π/4的值。

莱布尼茨级数的历史背景

莱布尼茨级数的发现是微积分学发展史上的一个重要里程碑。莱布尼茨与艾萨克·牛顿独立发现了微积分，他们的工作奠定了现代数学分析的基础。莱布尼茨级数不仅展示了无穷级数在计算π方面的应用，也反映了莱布尼茨对无穷小量和无穷级数的深刻理解。

莱布尼茨级数的应用

尽管莱布尼茨级数收敛速度慢，但它在数学教育中仍然有其价值。它提供了一个直观的方式来理解π的值，并展示了无穷级数在解决数学问题中的应用。此外，莱布尼茨级数也是研究π的数值方法的一个重要案例，它启发了后续更高效算法的发展。

结论

莱布尼茨级数是π的一个经典表示，它不仅在数学史上占有一席之地，而且在数学教育和研究中仍然具有启发性。

虽然它的收敛速度限制了它在实际计算中的应用，但它仍然是理解和探索π性质的一个有趣和有价值的工具。

通过莱布尼茨级数，我们可以窥见历史上数学家们对无穷概念的探索和应用，以及他们在数学分析领域的杰出贡献。

为什么莱布尼茨级数计算值接近 pi 呢？

简单点解释

想象你有一个完美的圆形，你想要知道这个圆的面积是多少。但是，你没有直接计算圆面积的工具，所以你决定用一种巧妙的方法来估计它。

你开始画一个正方形，这个正方形的边长和你圆的直径一样长。你知道正方形的面积是边长的平方，所以你很容易就算出了正方形的面积。但是，你实际上想要的是圆的面积，而不是正方形的面积。

这时，你注意到如果你把正方形的四个角切掉，剩下的形状会更接近圆形。

你开始想：“如果我切掉足够多的角，剩下的形状不就是圆了吗？”于是，你开始一个接一个地切掉正方形的角，每次都让剩下的形状更像圆。

每次切掉一个角，你都计算一下剩下的形状的面积。

你发现，每次切掉的角的面积都可以通过1除以一个奇数来表示（比如1/3、1/5、1/7等等）。你把这些面积加起来，然后从正方形的总面积里一个一个地减去这些值。

你发现，随着你切掉越来越多的角，剩下的形状的面积越来越接近圆的真实面积。你把这个不断逼近的过程变成了一个数学公式，就是莱布尼茨级数：

π/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...

这个级数告诉你，如果你继续这个过程，最终你会得到一个非常接近π/4的值。因为圆的面积是π乘以半径的平方，而半径在这里是1，所以这个级数实际上在帮你计算π。

所以，莱布尼茨级数通过不断地从正方形的面积中减去一些特定的值（这些值是1除以奇数），来逼近圆的面积。这就是为什么莱布尼茨级数能够计算出接近π的值的原因。

为啥一会加，一会儿减呢？

莱布尼茨级数中的加号和减号来源于对圆内接正多边形和外切正多边形面积的交替修正。

这个级数实际上是通过不断逼近的方式来计算圆的面积，从而得到π的值。

加法（正项）：当我们考虑圆内接正多边形时，随着边数的增加，多边形的边越来越接近圆的曲线，其面积也越来越接近圆的面积。但是，对于任何有限边数的多边形，它的面积总是小于圆的面积。因此，当我们计算内接多边形的面积并将其加入总和时，我们实际上是在逐步增加对圆面积的估计。
减法（负项）：另一方面，当我们考虑圆外切正多边形时，其面积总是大于圆的面积。为了修正这个过度估计，我们需要从当前的总和中减去外切多边形超出圆面积的部分。这就是为什么在级数中会出现负项，它们代表了从估计值中减去的部分。

莱布尼茨级数通过这种方式交替地增加和减少面积估计值，最终趋向于圆的真实面积。这个过程可以用以下的几何直观来理解：

我们从一个正方形开始，它的面积是边长的平方，即1^2 = 1。
我们切掉四个角，剩下的部分更接近圆形。切掉的每个角的面积可以用1除以一个奇数来表示，比如1/3、1/5等。这些切掉的部分是多边形面积的修正，所以我们用减法来表示。
然后，我们再添加一些较小的角回来，因为内接多边形的面积还是小于圆的面积。这些添加的部分也是用1除以奇数来表示，但是这次是奇数序列中的下一个数，比如1/7、1/9等。这些添加的部分是用加法来表示的。

通过这种交替的加减过程，我们逐渐逼近圆的面积，而整个级数的和就会趋向于π/4。

将这个和乘以4，我们就得到了π的估计值。这就是为什么莱布尼茨级数中会有加法和减法的原因。

官方的解释

莱布尼茨级数计算值接近π的原因在于它是基于圆的面积和周长的几何关系推导出来的。

莱布尼茨级数实际上是对圆的面积进行无穷级数展开的结果，这个级数的求和会收敛到π/4。

圆的面积公式是A = πr²，其中A是面积，r是半径。对于单位圆（半径为1），面积就是π。

另一方面，圆的周长公式是C = 2πr，对于单位圆，周长是2π。莱布尼茨级数利用了圆的周长和面积之间的关系来计算π。

莱布尼茨级数的几何解释可以从单位圆内接正多边形和外切正多边形的面积开始理解。随着正多边形边数的增加，内接正多边形的面积越来越接近圆的面积，而外切正多边形的面积则越来越大于圆的面积。当正多边形的边数趋于无穷大时，内接和外切正多边形的面积将分别收敛于圆的面积。

莱布尼茨级数通过交错相加和相减的方式来逼近这个极限。级数的每一项可以看作是计算内接和外切正多边形面积差异的一部分。当分母（即正多边形的边数）越来越大时，每一项的贡献越来越小，但它们加起来越来越接近圆的实际面积。

数学上，莱布尼茨级数是基于泰勒级数展开和交错级数的性质。级数中的每一项可以看作是函数f(x) = x^2在x = 0处的导数的泰勒级数展开的一部分。通过交错相加和相减，这个级数能够收敛到π/4。

因此，莱布尼茨级数之所以能够计算出一个接近π的值，是因为它基于圆的几何属性和无穷级数的数学原理。

尽管这个级数收敛得比较慢，但理论上，如果项数足够多，它的和将无限接近于π/4，从而可以用于计算π的近似值。

程序

我们用程序模拟一下:

package com.github.houbb.how.calc.pi;

// 作者：老马啸西风
public class LeibnizPiSeries {

    public static void main(String[] args) {
        int N = 1000000; // 级数的项数，越大越精确
        double pi = calculatePiUsingLeibnizSeries(N);
        System.out.printf("Approximation of Pi using Leibniz Series (N = %d): %.15f%n", N, pi);
    }

    public static double calculatePiUsingLeibnizSeries(int N) {
        double pi = 0.0;
        for (int i = 0; i < N; i++) {
            pi += Math.pow(-1, i) / (2 * i + 1);
        }
        pi *= 4; // 莱布尼茨级数计算的是π/4，所以要乘以4得到π
        return pi;
    }

}

我们模拟 100W 次，结果为：

Approximation of Pi using Leibniz Series (N = 1000000): 3.141591653589774

详细介绍一下楚德诺夫斯基算法（快速收敛算法）

楚德诺夫斯基算法是一种计算圆周率π的快速收敛算法，由美国数学家和工程师兄弟阿纳托利·楚德诺夫斯基（Anatoly Chudnovsky）和米哈伊尔·楚德诺夫斯基（Mikhail Chudnovsky）在20世纪60年代提出。这个算法以其极高的收敛速度和计算效率而著称，在计算π的极高位数时特别有效。

算法原理

楚德诺夫斯基算法基于一个快速收敛的超几何级数，该级数与拉马努金的某些公式有相似之处。算法的核心是一个关于π的表达式，它涉及到一个特殊的函数——拉马努金函数（Ramanujan’s function），以及一个复杂的数学恒等式。

算法的数学表达式

楚德诺夫斯基算法的数学表达式如下：

当然，楚德诺夫斯基算法的数学表达式可以用纯文本的形式表示如下：

1/pi = sum from k=0 to infinity of (64/log(2*pi))^(2k) * (2k)! / ((2k+1)(2k+2)(2k+3))

这个表达式定义了一个无穷级数，其中每一项都是一个分数，分子是双阶乘（(2k)!）乘以一个常数因子（(64/log(2*pi))^(2k)），分母是连续三个奇数的乘积（(2k+1)(2k+2)(2k+3)）。

通过计算这个级数的前几项并求和，我们可以得到π的近似值。随着项数的增加，这个和会越来越接近π的实际值。

这个级数的每一项都是一个分数，分母是关于k的函数，分子是一个固定的常数。随着k的增加，每一项的大小迅速减小，这意味着级数收敛得非常快。

算法的优势

楚德诺夫斯基算法的主要优势在于其收敛速度。相比于其他计算π的算法，如莱布尼茨级数或阿基米德方法，楚德诺夫斯基算法能够在计算较少项数的情况下获得非常高的精度。这使得它在计算π的极高位数时非常有用。

算法的应用

楚德诺夫斯基算法因其高效性，在计算圆周率的高精度数值时被广泛使用。它不仅在数学研究中有应用，还被用于各种需要高精度π值的科学和工程计算中。

计算记录

楚德诺夫斯基兄弟使用这个算法计算了超过十亿位的π值。

在计算机的帮助下，这个算法已经被用来计算π的数万亿位，创造了多项世界纪录。

价值

楚德诺夫斯基算法是计算圆周率π的一个非常有效的工具，它以其快速收敛的特性和高计算效率在数学和科学计算中占有重要地位。

通过这个算法，我们可以在有限的计算资源下，计算出π的极高精度数值，这对于科学研究和工程应用都有着重要的意义。

为什么这个公式接近于 pi 呢？

π的值与其在复平面上的解析性质密切相关。楚德诺夫斯基算法的公式揭示了π与某些特殊函数（如拉马努金函数）之间的联系。

这些函数在数学上有着深刻的内涵，它们在复分析、数论和其他数学分支中都有重要应用。

楚德诺夫斯基算法之所以能够有效地近似计算出π的值，并且具有快速收敛的特性，是因为它利用了解析几何和复分析中的深刻数学原理。

这个算法的基础是黎曼ζ函数（Riemann zeta function）和它的正则化形式，以及与π的某种内在联系。

黎曼ζ函数在复数域内有一个很好的性质，即它在复平面上的某些点（除了1这个极点之外）的值与π有关。

特别是，当ζ函数的参数s接近1时，它与π的关系变得尤为重要。通过适当的变换和简化，可以推导出一个与π相关的表达式，这就是楚德诺夫斯基算法中使用的级数形式。

计算公式：

楚德诺夫斯基算法基于解析表达式，特别是利用了圆周率的级数展开式。

算法的核心是一个快速收敛的级数，该级数的每一项都涉及到一个特定的数学函数——拉马努金函数（Ramanujan’s function）：

ζ(s)函数表达式可以表示为：

zeta(s) = 1/(s-1) + 7/8*(1/(s-1))^2 + 31/64*(1/(s-1))^3 + ...

这个表达式是拉马努金函数的一个特例，它是一个关于复变量s的幂级数，其中每一项都是1除以(s-1)的幂，幂的指数是递增的奇数。

这个级数在s接近1时收敛，并且在计算圆周率π时非常有用。

其中，zeta(s) 是黎曼ζ函数的正则化形式，它在计算圆周率时非常有用。

楚德诺夫斯基兄弟发现了一个更有效的公式，用于计算π的第n位数字：

pi = sum from k=0 to infinity of (-1)^k / (2k+1)^2 * (64/log(2*pi))^(2k)

这个表达式表示了一个无穷级数，其中每一项都是一个分数，分数的分子是(-1)的k次幂，分母是(2k+1)的平方，整个分数乘以(64/log(2*pi))^(2k)。

随着k的增加，级数中的每一项都会变得更小，使得级数快速收敛到π的值。

拉马努金函数

拉马努金函数是一个与π的计算有关的函数，它在某些π的级数展开中出现。

拉马努金函数与π的关系并不是直接通过楚德诺夫斯基算法体现的，但在其他类型的π的计算中，拉马努金函数提供了一种有效的方法来表示π的小数展开。

为什么收敛快？

楚德诺夫斯基算法中的级数收敛速度快的原因在于它的项具有特定的形式。

具体来说，每一项的形式是：

保证了项的符号交替变化。这意味着随着k的增加，每一项的贡献迅速减小，而且减小的速度非常快。因此，级数的前几项就包含了π值的大部分信息，而后续的项对总和的贡献微乎其微。

基本解释

楚德诺夫斯基算法之所以收敛速度快，主要归功于其级数项的特殊结构和数学性质。

以下是一些关键因素，解释了为什么这个算法能够快速收敛到π的值：

级数的结构：楚德诺夫斯基算法使用的级数是一个交错级数，其中的项具有交替的符号和特定的分母结构。每一项的形式是(-1)^k / (2k+1)^2，乘以一个随着k增加而指数级减小的因子(64/ log(2π))^(2k)。这种结构意味着随着k的增加，每一项的值迅速减小。
指数级减小的因子：级数中(64/ log(2π))^(2k)这一因子是关键，因为它使得每一项的系数以指数速率减小。这意味着级数的后期项对总和的贡献微乎其微，因此，只需要计算相对较少的项数就能达到非常高的精度。
数学优化：楚德诺夫斯基兄弟通过数学分析和优化，发现了这个级数，它在计算π时具有快速收敛的特性。这个级数的设计充分利用了数学中的某些恒等式和性质，从而实现了对π的高效计算。
收敛速度的理论基础：从数学理论上讲，一个级数收敛的速度与其项的衰减速率有关。楚德诺夫斯基算法中的级数项的衰减速率非常快，这是因为它们不仅仅是线性衰减或二次衰减，而是指数级的衰减。
实际计算效率：在实际计算中，楚德诺夫斯基算法的快速收敛特性意味着可以在较少的迭代次数内获得高精度的π值。这对于计算大量位数的π尤其重要，因为它可以显著减少所需的计算资源和时间。

综上所述，楚德诺夫斯基算法之所以收敛速度快，是因为它巧妙地利用了数学结构和性质，设计出了一个收敛速度极快的级数。这使得它成为了计算π值的一个非常有效的工具，尤其是在需要高精度计算时。

程序

package com.github.houbb.how.calc.pi;

import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.MathContext;
import java.math.RoundingMode;

// 作者：老马啸西风
// 参考：https://github.com/alanzheng88/JavaPI/blob/master/src/org/alanzheng/JavaPI.java
public class ChudnovskyPiApproximation {

    static final BigInteger BIG_A = BigInteger.valueOf(13591409);
    static final BigInteger BIG_B = BigInteger.valueOf(545140134);
    static final BigInteger BIG_C = new BigInteger("-262537412640768000");
    static final BigInteger THREE = BigInteger.valueOf(3);
    static final BigInteger SIX = BigInteger.valueOf(6);
    static final BigInteger TWELVE = BigInteger.valueOf(12);
    static final BigInteger SIXTEEN = BigInteger.valueOf(16);

    // complexity of algorithm O(n * log(n)^3)
    public static BigDecimal getValueOfPi(int num) {
        BigDecimal C = new BigDecimal(426880 * Math.sqrt(10005));
        MathContext mc = new MathContext(num+1);
        BigDecimal sum = new BigDecimal(BIG_A);
        BigInteger BIG_D;
        BigInteger L = BIG_A;
        BigInteger X = BigInteger.ONE;
        BigDecimal M = BigDecimal.ONE;
        BigInteger K = SIX;
        for (int d = 1; d <= num; d++) {
            BIG_D = BigInteger.valueOf(Long.valueOf(d));
            L = L.add(BIG_B);
            X = X.multiply(BIG_C);
            M = M.multiply( new BigDecimal(K.pow(3).subtract(SIXTEEN.multiply(K)))
                    .divide(new BigDecimal(BIG_D)
                            .add(BigDecimal.ONE).pow(3), mc)  );
            K = K.add(TWELVE);
            sum = sum.add( M.multiply(new BigDecimal(L))
                    .divide(new BigDecimal(X), mc) );
        }
        BigDecimal PI = C.divide(sum, mc);
        return PI.setScale(num, RoundingMode.FLOOR);
    }

    public static void main(String[] args) {
        System.out.println(getValueOfPi(100));
    }

}

结果：

3.1415926535897413249367471616497999120156209473408921170808260143385335751657450817201892302297784403

速度确实非常快，不过这个结果是存疑的。

算法随便从网上找了一份，且存在精度丢失等问题，计算的不见得准确。

比如网上有前 100 位，可以简单对比下：

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679

参考资料

前言
历史与发展
- 回头看中国
- 抬头看世界
阿基米德方法（几何方法）
- 阿基米德的穷竭法
蒙特卡洛方法（统计方法）
- 说明
- 程序代码
莱布尼茨级数
详细介绍一下楚德诺夫斯基算法（快速收敛算法）
更多 pi 解法
- 拓展
- 小结
参考资料

前言

历史与发展

回头看中国

早期记载与探索

祖冲之的贡献

后世的影响与纪念

抬头看世界

1. 阿基米德方法（几何方法）

2. 蒙特卡洛方法（统计方法）

3. 莱布尼茨级数（解析方法）

4. 楚德诺夫斯基算法（快速收敛算法）

阿基米德方法（几何方法）

阿基米德的穷竭法

蒙特卡洛方法（统计方法）

说明

蒙特卡洛方法的基本原理

蒙特卡洛方法的实现

蒙特卡洛方法的优点

蒙特卡洛方法的局限性

程序代码

莱布尼茨级数

说明

莱布尼茨级数的收敛性

莱布尼茨级数的历史背景

莱布尼茨级数的应用

结论

为什么莱布尼茨级数计算值接近 pi 呢？

简单点解释

为啥一会加，一会儿减呢？

官方的解释

程序

详细介绍一下 楚德诺夫斯基算法（快速收敛算法）

算法原理

算法的数学表达式

算法的优势

算法的应用

计算记录

价值

为什么这个公式接近于 pi 呢？

拉马努金函数

为什么收敛快？

基本解释

程序

更多 pi 解法

拓展

高斯积分（Gaussian Integral）

布冯针问题（Buffon’s Needle）

贝塞尔问题（Bessel’s Problem）

沃利斯公式（Wallis’s Formula）

马青公式（Machin’s formula）

格雷戈里-莱布尼茨级数（Gregory-Leibniz series）

正弦级数

切比雪夫公式

芝诺级数

拉马努金公式

快速傅里叶变换（FFT）

迭代方法

数值积分方法

小结

参考资料

更多学习

详细介绍一下楚德诺夫斯基算法（快速收敛算法）