我正在寻找最快的方法来获得π的值,作为一个个人挑战。更具体地说,我使用的方法不涉及使用#define常量M_PI,或硬编码的数字。

下面的程序测试了我所知道的各种方法。从理论上讲,内联汇编版本是最快的选择,尽管显然不能移植。我将它作为一个基准,与其他版本进行比较。在我的测试中,使用内置函数,4 * atan(1)版本在GCC 4.2上是最快的,因为它自动将atan(1)折叠成一个常量。通过指定-fno-builtin, atan2(0, -1)版本是最快的。

下面是主要的测试程序(pitimes.c):

#include <math.h>
#include <stdio.h>
#include <time.h>

#define ITERS 10000000
#define TESTWITH(x) {                                                       \
    diff = 0.0;                                                             \
    time1 = clock();                                                        \
    for (i = 0; i < ITERS; ++i)                                             \
        diff += (x) - M_PI;                                                 \
    time2 = clock();                                                        \
    printf("%s\t=> %e, time => %f\n", #x, diff, diffclock(time2, time1));   \
}

static inline double
diffclock(clock_t time1, clock_t time0)
{
    return (double) (time1 - time0) / CLOCKS_PER_SEC;
}

int
main()
{
    int i;
    clock_t time1, time2;
    double diff;

    /* Warmup. The atan2 case catches GCC's atan folding (which would
     * optimise the ``4 * atan(1) - M_PI'' to a no-op), if -fno-builtin
     * is not used. */
    TESTWITH(4 * atan(1))
    TESTWITH(4 * atan2(1, 1))

#if defined(__GNUC__) && (defined(__i386__) || defined(__amd64__))
    extern double fldpi();
    TESTWITH(fldpi())
#endif

    /* Actual tests start here. */
    TESTWITH(atan2(0, -1))
    TESTWITH(acos(-1))
    TESTWITH(2 * asin(1))
    TESTWITH(4 * atan2(1, 1))
    TESTWITH(4 * atan(1))

    return 0;
}

内联汇编的东西(fldpi.c)只适用于x86和x64系统:

double
fldpi()
{
    double pi;
    asm("fldpi" : "=t" (pi));
    return pi;
}

和一个构建脚本,构建我正在测试的所有配置(build.sh):

#!/bin/sh
gcc -O3 -Wall -c           -m32 -o fldpi-32.o fldpi.c
gcc -O3 -Wall -c           -m64 -o fldpi-64.o fldpi.c

gcc -O3 -Wall -ffast-math  -m32 -o pitimes1-32 pitimes.c fldpi-32.o
gcc -O3 -Wall              -m32 -o pitimes2-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -fno-builtin -m32 -o pitimes3-32 pitimes.c fldpi-32.o -lm
gcc -O3 -Wall -ffast-math  -m64 -o pitimes1-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall              -m64 -o pitimes2-64 pitimes.c fldpi-64.o -lm
gcc -O3 -Wall -fno-builtin -m64 -o pitimes3-64 pitimes.c fldpi-64.o -lm

除了在各种编译器标志之间进行测试(我也比较了32位和64位,因为优化是不同的),我还尝试切换测试的顺序。但是,atan2(0, -1)版本在每次测试中仍然名列前茅。


当前回答

Chudnovsky算法非常快如果你不介意做一个平方根和几个逆运算的话。它在2次迭代中收敛到两倍精度。

/*
    Chudnovsky algorithm for computing PI
*/

#include <iostream>
#include <cmath>
using namespace std;

double calc_PI(int K=2) {

    static const int A = 545140134;
    static const int B = 13591409;
    static const int D = 640320;

    const double ID3 = 1./ (double(D)*double(D)*double(D));

    double sum = 0.;
    double b   = sqrt(ID3);
    long long int p = 1;
    long long int a = B;

    sum += double(p) * double(a)* b;

    // 2 iterations enough for double convergence
    for (int k=1; k<K; ++k) {
        // A*k + B
        a += A;
        // update denominator
        b *= ID3;
        // p = (-1)^k 6k! / 3k! k!^3
        p *= (6*k)*(6*k-1)*(6*k-2)*(6*k-3)*(6*k-4)*(6*k-5);
        p /= (3*k)*(3*k-1)*(3*k-2) * k*k*k;
        p = -p;

        sum += double(p) * double(a)* b;
    }

    return 1./(12*sum);
}

int main() {

    cout.precision(16);
    cout.setf(ios::fixed);

    for (int k=1; k<=5; ++k) cout << "k = " << k << "   PI = " << calc_PI(k) << endl;

    return 0;
}

结果:

k = 1   PI = 3.1415926535897341
k = 2   PI = 3.1415926535897931
k = 3   PI = 3.1415926535897931
k = 4   PI = 3.1415926535897931
k = 5   PI = 3.1415926535897931

其他回答

如果你想计算π值的近似值(出于某种原因),你应该尝试二进制提取算法。Bellard对BBP的改进给出了O(N²)中的PI。


如果你想获得π值的近似值来进行计算,那么:

PI = 3.141592654

当然,这只是一个近似值,并不完全准确。误差略大于0.00000000004102。(4个十万亿分之一,大约4/10,000,000,000)。


如果你想用π做数学运算,那就准备好铅笔和纸,或者电脑代数包,然后使用π的精确值π。

如果你真的想要一个公式,这个很有趣:

π = -i ln(-1)

下面是我在高中时学过的计算圆周率的技巧。

我之所以分享它,是因为我认为它足够简单,任何人都可以无限期地记住它,而且它教会了你“蒙特卡罗”方法的概念——这是一种统计方法,可以得到答案,这些答案不会立即通过随机过程演绎出来。

画一个正方形,在这个正方形内画一个象限(半圆的四分之一)(一个半径等于正方形边的象限,这样它就能尽可能多地填充正方形)

现在向正方形投掷飞镖,并记录飞镖落在何处——也就是说,在正方形内任意选择一个点。当然,它落在了正方形内部,但它落在半圆内部吗?记录这个事实。

重复此过程多次,你会发现半圆内的点数量与抛出的总数量之比为x。

由于正方形的面积是r乘以r,可以推导出半圆的面积是x乘以r乘以r(即x乘以r的平方)。因此x乘以4会得到。

这不是一个快速使用的方法。但这是蒙特卡罗方法的一个很好的例子。如果你环顾四周,你可能会发现许多超出你计算能力的问题都可以用这种方法来解决。

蒙特卡罗方法,如前所述,应用了一些伟大的概念,但很明显,它不是最快的,不是从任何合理的标准来看。此外,这完全取决于你想要什么样的准确性。我所知道的最快的π是数字硬编码的π。看看圆周率和圆周率,有很多公式。

Here is a method that converges quickly — about 14 digits per iteration. PiFast, the current fastest application, uses this formula with the FFT. I'll just write the formula, since the code is straightforward. This formula was almost found by Ramanujan and discovered by Chudnovsky. It is actually how he calculated several billion digits of the number — so it isn't a method to disregard. The formula will overflow quickly and, since we are dividing factorials, it would be advantageous then to delay such calculations to remove terms.

在那里,

下面是Brent-Salamin算法。维基百科提到,当a和b“足够接近”时,(a + b)²/ 4t将是π的近似值。我不确定“足够接近”是什么意思,但从我的测试来看,一次迭代得到2位数字,两次得到7位,3次得到15位,当然这是双精度,所以它可能会有一个基于它的表示的错误,真实的计算可能会更准确。

let pi_2 iters =
    let rec loop_ a b t p i =
        if i = 0 then a,b,t,p
        else
            let a_n = (a +. b) /. 2.0 
            and b_n = sqrt (a*.b)
            and p_n = 2.0 *. p in
            let t_n = t -. (p *. (a -. a_n) *. (a -. a_n)) in
            loop_ a_n b_n t_n p_n (i - 1)
    in 
    let a,b,t,p = loop_ (1.0) (1.0 /. (sqrt 2.0)) (1.0/.4.0) (1.0) iters in
    (a +. b) *. (a +. b) /. (4.0 *. t)

最后,来点圆周率高尔夫(800位数字)怎么样?160个字符!

int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

比GMPY2和MPmath内置更快:45分钟十亿:


我尝试了几种方法;Manchin, AGM和Chudnovsky兄弟。Chudnovsky和Binary Split是最快的: 我的github: https://github.com/Overboard-code/Pi-Pourri

我的Binary Split Chudnovsky的速度大约是内置gmpy2.const_pi()的两倍。MPmath.mp.pi()计算10亿需要50分钟,所以它几乎和Chudnovsky一样快。

我也非常感谢表演技巧。我不确定我的代码是否完美。它是100%准确的(所有公式都同意1亿),但也许可以更快?

我尝试了gmpy2.const_pi()到1亿个数字,在同一台机器上,Chudnovsky花了300秒,而Chudnovsky花了150秒。Pi.txt和pi2.txt是一样的。

在不到一个小时的时间里,我在我的旧i7 16GB笔记本电脑上输入了10亿个数字。

以下是我尝试过的12种方法中最快的一种:

class PiChudnovsky:
    """Version of Chudnovsky Bros using Binary Splitting 
        So far this is the winner for fastest time to a million digits on my older intel i7
    """
    A = mpz(13591409)
    B = mpz(545140134)
    C = mpz(640320)
    D = mpz(426880)
    E = mpz(10005)
    C3_24  = pow(C, mpz(3)) // mpz(24)
    #DIGITS_PER_TERM = math.log(53360 ** 3) / math.log(10)  #=> 14.181647462725476
    DIGITS_PER_TERM = 14.181647462725476
    MMILL = mpz(1000000)

    def __init__(self,ndigits):
        """ Initialization
        :param int ndigits: digits of PI computation
        """
        self.ndigits = ndigits
        self.n      = mpz(self.ndigits // self.DIGITS_PER_TERM + 1)
        self.prec   = mpz((self.ndigits + 1) * LOG2_10)
        self.one_sq = pow(mpz(10),mpz(2 * ndigits))
        self.sqrt_c = isqrt(self.E * self.one_sq)
        self.iters  = mpz(0)
        self.start_time = 0

    def compute(self):
        """ Computation """
        try:
            self.start_time = time.time()
            logging.debug("Starting {} formula to {:,} decimal places"
                .format(name,ndigits) )
            __, q, t = self.__bs(mpz(0), self.n)  # p is just for recursion
            pi = (q * self.D * self.sqrt_c) // t
            logging.debug('{} calulation Done! {:,} iterations and {:.2f} seconds.'
                .format( name, int(self.iters),time.time() - self.start_time))
            get_context().precision= int((self.ndigits+10) * LOG2_10)
            pi_s = pi.digits() # digits() gmpy2 creates a string 
            pi_o = pi_s[:1] + "." + pi_s[1:]
            return pi_o,int(self.iters),time.time() - self.start_time
        except Exception as e:
            print (e.message, e.args)
            raise

    def __bs(self, a, b):
        """ PQT computation by BSA(= Binary Splitting Algorithm)
        :param int a: positive integer
        :param int b: positive integer
        :return list [int p_ab, int q_ab, int t_ab]
        """
        try:
            self.iters += mpz(1)
            if self.iters % self.MMILL  == mpz(0):
                logging.debug('Chudnovsky ... {:,} iterations and {:.2f} seconds.'
                    .format( int(self.iters),time.time() - self.start_time))
            if a + mpz(1) == b:
                if a == mpz(0):
                    p_ab = q_ab = mpz(1)
                else:
                    p_ab = mpz((mpz(6) * a - mpz(5)) * (mpz(2) * a - mpz(1)) * (mpz(6) * a - mpz(1)))
                    q_ab = pow(a,mpz(3)) * self.C3_24
                t_ab = p_ab * (self.A + self.B * a)
                if a & 1:
                    t_ab *= mpz(-1)
            else:
                m = (a + b) // mpz(2)
                p_am, q_am, t_am = self.__bs(a, m)
                p_mb, q_mb, t_mb = self.__bs(m, b)
                p_ab = p_am * p_mb
                q_ab = q_am * q_mb
                t_ab = q_mb * t_am + p_am * t_mb
            return [p_ab, q_ab, t_ab]
        except Exception as e:
            print (e.message, e.args)
            raise

以下是在45分钟内输出的10亿位数:

python pi-pourri.py -v -d 1,000,000,000 -a 10 

[INFO] 2022-10-03 09:22:51,860 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-03 09:25:00,543 compute: MainProcess Starting   Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-03 09:25:04,995 __bs: MainProcess Chudnovsky ... 1,000,000 iterations and 4.45 seconds.
[DEBUG] 2022-10-03 09:25:10,836 __bs: MainProcess Chudnovsky ... 2,000,000 iterations and 10.29 seconds.
[DEBUG] 2022-10-03 09:25:18,227 __bs: MainProcess Chudnovsky ... 3,000,000 iterations and 17.68 seconds.
[DEBUG] 2022-10-03 09:25:24,512 __bs: MainProcess Chudnovsky ... 4,000,000 iterations and 23.97 seconds.
[DEBUG] 2022-10-03 09:25:35,670 __bs: MainProcess Chudnovsky ... 5,000,000 iterations and 35.13 seconds.
[DEBUG] 2022-10-03 09:25:41,376 __bs: MainProcess Chudnovsky ... 6,000,000 iterations and 40.83 seconds.
[DEBUG] 2022-10-03 09:25:49,238 __bs: MainProcess Chudnovsky ... 7,000,000 iterations and 48.69 seconds.
[DEBUG] 2022-10-03 09:25:55,646 __bs: MainProcess Chudnovsky ... 8,000,000 iterations and 55.10 seconds.
[DEBUG] 2022-10-03 09:26:15,043 __bs: MainProcess Chudnovsky ... 9,000,000 iterations and 74.50 seconds.
[DEBUG] 2022-10-03 09:26:21,437 __bs: MainProcess Chudnovsky ... 10,000,000 iterations and 80.89 seconds.
[DEBUG] 2022-10-03 09:26:26,587 __bs: MainProcess Chudnovsky ... 11,000,000 iterations and 86.04 seconds.
[DEBUG] 2022-10-03 09:26:34,777 __bs: MainProcess Chudnovsky ... 12,000,000 iterations and 94.23 seconds.
[DEBUG] 2022-10-03 09:26:41,231 __bs: MainProcess Chudnovsky ... 13,000,000 iterations and 100.69 seconds.
[DEBUG] 2022-10-03 09:26:52,972 __bs: MainProcess Chudnovsky ... 14,000,000 iterations and 112.43 seconds.
[DEBUG] 2022-10-03 09:26:59,517 __bs: MainProcess Chudnovsky ... 15,000,000 iterations and 118.97 seconds.
[DEBUG] 2022-10-03 09:27:07,932 __bs: MainProcess Chudnovsky ... 16,000,000 iterations and 127.39 seconds.
[DEBUG] 2022-10-03 09:27:14,036 __bs: MainProcess Chudnovsky ... 17,000,000 iterations and 133.49 seconds.
[DEBUG] 2022-10-03 09:27:51,629 __bs: MainProcess Chudnovsky ... 18,000,000 iterations and 171.09 seconds.
[DEBUG] 2022-10-03 09:27:58,176 __bs: MainProcess Chudnovsky ... 19,000,000 iterations and 177.63 seconds.
[DEBUG] 2022-10-03 09:28:06,704 __bs: MainProcess Chudnovsky ... 20,000,000 iterations and 186.16 seconds.
[DEBUG] 2022-10-03 09:28:13,376 __bs: MainProcess Chudnovsky ... 21,000,000 iterations and 192.83 seconds.
[DEBUG] 2022-10-03 09:28:18,737 __bs: MainProcess Chudnovsky ... 22,000,000 iterations and 198.19 seconds.
[DEBUG] 2022-10-03 09:28:31,095 __bs: MainProcess Chudnovsky ... 23,000,000 iterations and 210.55 seconds.
[DEBUG] 2022-10-03 09:28:37,789 __bs: MainProcess Chudnovsky ... 24,000,000 iterations and 217.25 seconds.
[DEBUG] 2022-10-03 09:28:46,171 __bs: MainProcess Chudnovsky ... 25,000,000 iterations and 225.63 seconds.
[DEBUG] 2022-10-03 09:28:52,933 __bs: MainProcess Chudnovsky ... 26,000,000 iterations and 232.39 seconds.
[DEBUG] 2022-10-03 09:29:13,524 __bs: MainProcess Chudnovsky ... 27,000,000 iterations and 252.98 seconds.
[DEBUG] 2022-10-03 09:29:19,676 __bs: MainProcess Chudnovsky ... 28,000,000 iterations and 259.13 seconds.
[DEBUG] 2022-10-03 09:29:28,196 __bs: MainProcess Chudnovsky ... 29,000,000 iterations and 267.65 seconds.
[DEBUG] 2022-10-03 09:29:34,720 __bs: MainProcess Chudnovsky ... 30,000,000 iterations and 274.18 seconds.
[DEBUG] 2022-10-03 09:29:47,075 __bs: MainProcess Chudnovsky ... 31,000,000 iterations and 286.53 seconds.
[DEBUG] 2022-10-03 09:29:53,746 __bs: MainProcess Chudnovsky ... 32,000,000 iterations and 293.20 seconds.
[DEBUG] 2022-10-03 09:29:59,099 __bs: MainProcess Chudnovsky ... 33,000,000 iterations and 298.56 seconds.
[DEBUG] 2022-10-03 09:30:07,511 __bs: MainProcess Chudnovsky ... 34,000,000 iterations and 306.97 seconds.
[DEBUG] 2022-10-03 09:30:14,279 __bs: MainProcess Chudnovsky ... 35,000,000 iterations and 313.74 seconds.
[DEBUG] 2022-10-03 09:31:31,710 __bs: MainProcess Chudnovsky ... 36,000,000 iterations and 391.17 seconds.
[DEBUG] 2022-10-03 09:31:38,454 __bs: MainProcess Chudnovsky ... 37,000,000 iterations and 397.91 seconds.
[DEBUG] 2022-10-03 09:31:46,437 __bs: MainProcess Chudnovsky ... 38,000,000 iterations and 405.89 seconds.
[DEBUG] 2022-10-03 09:31:53,285 __bs: MainProcess Chudnovsky ... 39,000,000 iterations and 412.74 seconds.
[DEBUG] 2022-10-03 09:32:05,602 __bs: MainProcess Chudnovsky ... 40,000,000 iterations and 425.06 seconds.
[DEBUG] 2022-10-03 09:32:12,220 __bs: MainProcess Chudnovsky ... 41,000,000 iterations and 431.68 seconds.
[DEBUG] 2022-10-03 09:32:20,708 __bs: MainProcess Chudnovsky ... 42,000,000 iterations and 440.17 seconds.
[DEBUG] 2022-10-03 09:32:27,552 __bs: MainProcess Chudnovsky ... 43,000,000 iterations and 447.01 seconds.
[DEBUG] 2022-10-03 09:32:32,986 __bs: MainProcess Chudnovsky ... 44,000,000 iterations and 452.44 seconds.
[DEBUG] 2022-10-03 09:32:53,904 __bs: MainProcess Chudnovsky ... 45,000,000 iterations and 473.36 seconds.
[DEBUG] 2022-10-03 09:33:00,832 __bs: MainProcess Chudnovsky ... 46,000,000 iterations and 480.29 seconds.
[DEBUG] 2022-10-03 09:33:09,198 __bs: MainProcess Chudnovsky ... 47,000,000 iterations and 488.66 seconds.
[DEBUG] 2022-10-03 09:33:16,000 __bs: MainProcess Chudnovsky ... 48,000,000 iterations and 495.46 seconds.
[DEBUG] 2022-10-03 09:33:27,921 __bs: MainProcess Chudnovsky ... 49,000,000 iterations and 507.38 seconds.
[DEBUG] 2022-10-03 09:33:34,778 __bs: MainProcess Chudnovsky ... 50,000,000 iterations and 514.24 seconds.
[DEBUG] 2022-10-03 09:33:43,298 __bs: MainProcess Chudnovsky ... 51,000,000 iterations and 522.76 seconds.
[DEBUG] 2022-10-03 09:33:49,959 __bs: MainProcess Chudnovsky ... 52,000,000 iterations and 529.42 seconds.
[DEBUG] 2022-10-03 09:34:29,294 __bs: MainProcess Chudnovsky ... 53,000,000 iterations and 568.75 seconds.
[DEBUG] 2022-10-03 09:34:36,176 __bs: MainProcess Chudnovsky ... 54,000,000 iterations and 575.63 seconds.
[DEBUG] 2022-10-03 09:34:41,576 __bs: MainProcess Chudnovsky ... 55,000,000 iterations and 581.03 seconds.
[DEBUG] 2022-10-03 09:34:50,161 __bs: MainProcess Chudnovsky ... 56,000,000 iterations and 589.62 seconds.
[DEBUG] 2022-10-03 09:34:56,811 __bs: MainProcess Chudnovsky ... 57,000,000 iterations and 596.27 seconds.
[DEBUG] 2022-10-03 09:35:09,382 __bs: MainProcess Chudnovsky ... 58,000,000 iterations and 608.84 seconds.
[DEBUG] 2022-10-03 09:35:16,206 __bs: MainProcess Chudnovsky ... 59,000,000 iterations and 615.66 seconds.
[DEBUG] 2022-10-03 09:35:24,295 __bs: MainProcess Chudnovsky ... 60,000,000 iterations and 623.75 seconds.
[DEBUG] 2022-10-03 09:35:31,095 __bs: MainProcess Chudnovsky ... 61,000,000 iterations and 630.55 seconds.
[DEBUG] 2022-10-03 09:35:52,139 __bs: MainProcess Chudnovsky ... 62,000,000 iterations and 651.60 seconds.
[DEBUG] 2022-10-03 09:35:58,781 __bs: MainProcess Chudnovsky ... 63,000,000 iterations and 658.24 seconds.
[DEBUG] 2022-10-03 09:36:07,399 __bs: MainProcess Chudnovsky ... 64,000,000 iterations and 666.86 seconds.
[DEBUG] 2022-10-03 09:36:12,847 __bs: MainProcess Chudnovsky ... 65,000,000 iterations and 672.30 seconds.
[DEBUG] 2022-10-03 09:36:19,763 __bs: MainProcess Chudnovsky ... 66,000,000 iterations and 679.22 seconds.
[DEBUG] 2022-10-03 09:36:32,351 __bs: MainProcess Chudnovsky ... 67,000,000 iterations and 691.81 seconds.
[DEBUG] 2022-10-03 09:36:39,078 __bs: MainProcess Chudnovsky ... 68,000,000 iterations and 698.53 seconds.
[DEBUG] 2022-10-03 09:36:47,830 __bs: MainProcess Chudnovsky ... 69,000,000 iterations and 707.29 seconds.
[DEBUG] 2022-10-03 09:36:54,701 __bs: MainProcess Chudnovsky ... 70,000,000 iterations and 714.16 seconds.
[DEBUG] 2022-10-03 09:39:39,357 __bs: MainProcess Chudnovsky ... 71,000,000 iterations and 878.81 seconds.
[DEBUG] 2022-10-03 09:39:46,199 __bs: MainProcess Chudnovsky ... 72,000,000 iterations and 885.66 seconds.
[DEBUG] 2022-10-03 09:39:54,956 __bs: MainProcess Chudnovsky ... 73,000,000 iterations and 894.41 seconds.
[DEBUG] 2022-10-03 09:40:01,639 __bs: MainProcess Chudnovsky ... 74,000,000 iterations and 901.10 seconds.
[DEBUG] 2022-10-03 09:40:14,219 __bs: MainProcess Chudnovsky ... 75,000,000 iterations and 913.68 seconds.
[DEBUG] 2022-10-03 09:40:19,680 __bs: MainProcess Chudnovsky ... 76,000,000 iterations and 919.14 seconds.
[DEBUG] 2022-10-03 09:40:26,625 __bs: MainProcess Chudnovsky ... 77,000,000 iterations and 926.08 seconds.
[DEBUG] 2022-10-03 09:40:35,212 __bs: MainProcess Chudnovsky ... 78,000,000 iterations and 934.67 seconds.
[DEBUG] 2022-10-03 09:40:41,914 __bs: MainProcess Chudnovsky ... 79,000,000 iterations and 941.37 seconds.
[DEBUG] 2022-10-03 09:41:03,218 __bs: MainProcess Chudnovsky ... 80,000,000 iterations and 962.68 seconds.
[DEBUG] 2022-10-03 09:41:10,213 __bs: MainProcess Chudnovsky ... 81,000,000 iterations and 969.67 seconds.
[DEBUG] 2022-10-03 09:41:18,344 __bs: MainProcess Chudnovsky ... 82,000,000 iterations and 977.80 seconds.
[DEBUG] 2022-10-03 09:41:25,261 __bs: MainProcess Chudnovsky ... 83,000,000 iterations and 984.72 seconds.
[DEBUG] 2022-10-03 09:41:37,663 __bs: MainProcess Chudnovsky ... 84,000,000 iterations and 997.12 seconds.
[DEBUG] 2022-10-03 09:41:44,680 __bs: MainProcess Chudnovsky ... 85,000,000 iterations and 1004.14 seconds.
[DEBUG] 2022-10-03 09:41:53,411 __bs: MainProcess Chudnovsky ... 86,000,000 iterations and 1012.87 seconds.
[DEBUG] 2022-10-03 09:41:58,926 __bs: MainProcess Chudnovsky ... 87,000,000 iterations and 1018.38 seconds.
[DEBUG] 2022-10-03 09:42:05,858 __bs: MainProcess Chudnovsky ... 88,000,000 iterations and 1025.32 seconds.
[DEBUG] 2022-10-03 09:42:46,163 __bs: MainProcess Chudnovsky ... 89,000,000 iterations and 1065.62 seconds.
[DEBUG] 2022-10-03 09:42:53,054 __bs: MainProcess Chudnovsky ... 90,000,000 iterations and 1072.51 seconds.
[DEBUG] 2022-10-03 09:43:02,030 __bs: MainProcess Chudnovsky ... 91,000,000 iterations and 1081.49 seconds.
[DEBUG] 2022-10-03 09:43:09,192 __bs: MainProcess Chudnovsky ... 92,000,000 iterations and 1088.65 seconds.
[DEBUG] 2022-10-03 09:43:21,533 __bs: MainProcess Chudnovsky ... 93,000,000 iterations and 1100.99 seconds.
[DEBUG] 2022-10-03 09:43:28,643 __bs: MainProcess Chudnovsky ... 94,000,000 iterations and 1108.10 seconds.
[DEBUG] 2022-10-03 09:43:37,372 __bs: MainProcess Chudnovsky ... 95,000,000 iterations and 1116.83 seconds.
[DEBUG] 2022-10-03 09:43:44,558 __bs: MainProcess Chudnovsky ... 96,000,000 iterations and 1124.02 seconds.
[DEBUG] 2022-10-03 09:44:06,555 __bs: MainProcess Chudnovsky ... 97,000,000 iterations and 1146.01 seconds.
[DEBUG] 2022-10-03 09:44:12,220 __bs: MainProcess Chudnovsky ... 98,000,000 iterations and 1151.68 seconds.
[DEBUG] 2022-10-03 09:44:19,278 __bs: MainProcess Chudnovsky ... 99,000,000 iterations and 1158.74 seconds.
[DEBUG] 2022-10-03 09:44:28,323 __bs: MainProcess Chudnovsky ... 100,000,000 iterations and 1167.78 seconds.
[DEBUG] 2022-10-03 09:44:35,211 __bs: MainProcess Chudnovsky ... 101,000,000 iterations and 1174.67 seconds.
[DEBUG] 2022-10-03 09:44:48,331 __bs: MainProcess Chudnovsky ... 102,000,000 iterations and 1187.79 seconds.
[DEBUG] 2022-10-03 09:44:54,835 __bs: MainProcess Chudnovsky ... 103,000,000 iterations and 1194.29 seconds.
[DEBUG] 2022-10-03 09:45:03,869 __bs: MainProcess Chudnovsky ... 104,000,000 iterations and 1203.33 seconds.
[DEBUG] 2022-10-03 09:45:10,967 __bs: MainProcess Chudnovsky ... 105,000,000 iterations and 1210.42 seconds.
[DEBUG] 2022-10-03 09:46:32,760 __bs: MainProcess Chudnovsky ... 106,000,000 iterations and 1292.22 seconds.
[DEBUG] 2022-10-03 09:46:39,872 __bs: MainProcess Chudnovsky ... 107,000,000 iterations and 1299.33 seconds.
[DEBUG] 2022-10-03 09:46:48,948 __bs: MainProcess Chudnovsky ... 108,000,000 iterations and 1308.41 seconds.
[DEBUG] 2022-10-03 09:46:54,611 __bs: MainProcess Chudnovsky ... 109,000,000 iterations and 1314.07 seconds.
[DEBUG] 2022-10-03 09:47:01,727 __bs: MainProcess Chudnovsky ... 110,000,000 iterations and 1321.18 seconds.
[DEBUG] 2022-10-03 09:47:14,525 __bs: MainProcess Chudnovsky ... 111,000,000 iterations and 1333.98 seconds.
[DEBUG] 2022-10-03 09:47:21,682 __bs: MainProcess Chudnovsky ... 112,000,000 iterations and 1341.14 seconds.
[DEBUG] 2022-10-03 09:47:30,610 __bs: MainProcess Chudnovsky ... 113,000,000 iterations and 1350.07 seconds.
[DEBUG] 2022-10-03 09:47:37,176 __bs: MainProcess Chudnovsky ... 114,000,000 iterations and 1356.63 seconds.
[DEBUG] 2022-10-03 09:47:59,642 __bs: MainProcess Chudnovsky ... 115,000,000 iterations and 1379.10 seconds.
[DEBUG] 2022-10-03 09:48:06,702 __bs: MainProcess Chudnovsky ... 116,000,000 iterations and 1386.16 seconds.
[DEBUG] 2022-10-03 09:48:15,483 __bs: MainProcess Chudnovsky ... 117,000,000 iterations and 1394.94 seconds.
[DEBUG] 2022-10-03 09:48:22,537 __bs: MainProcess Chudnovsky ... 118,000,000 iterations and 1401.99 seconds.
[DEBUG] 2022-10-03 09:48:35,714 __bs: MainProcess Chudnovsky ... 119,000,000 iterations and 1415.17 seconds.
[DEBUG] 2022-10-03 09:48:41,321 __bs: MainProcess Chudnovsky ... 120,000,000 iterations and 1420.78 seconds.
[DEBUG] 2022-10-03 09:48:48,408 __bs: MainProcess Chudnovsky ... 121,000,000 iterations and 1427.87 seconds.
[DEBUG] 2022-10-03 09:48:57,138 __bs: MainProcess Chudnovsky ... 122,000,000 iterations and 1436.60 seconds.
[DEBUG] 2022-10-03 09:49:04,328 __bs: MainProcess Chudnovsky ... 123,000,000 iterations and 1443.79 seconds.
[DEBUG] 2022-10-03 09:49:46,274 __bs: MainProcess Chudnovsky ... 124,000,000 iterations and 1485.73 seconds.
[DEBUG] 2022-10-03 09:49:52,833 __bs: MainProcess Chudnovsky ... 125,000,000 iterations and 1492.29 seconds.
[DEBUG] 2022-10-03 09:50:01,786 __bs: MainProcess Chudnovsky ... 126,000,000 iterations and 1501.24 seconds.
[DEBUG] 2022-10-03 09:50:08,975 __bs: MainProcess Chudnovsky ... 127,000,000 iterations and 1508.43 seconds.
[DEBUG] 2022-10-03 09:50:21,850 __bs: MainProcess Chudnovsky ... 128,000,000 iterations and 1521.31 seconds.
[DEBUG] 2022-10-03 09:50:28,962 __bs: MainProcess Chudnovsky ... 129,000,000 iterations and 1528.42 seconds.
[DEBUG] 2022-10-03 09:50:34,594 __bs: MainProcess Chudnovsky ... 130,000,000 iterations and 1534.05 seconds.
[DEBUG] 2022-10-03 09:50:43,647 __bs: MainProcess Chudnovsky ... 131,000,000 iterations and 1543.10 seconds.
[DEBUG] 2022-10-03 09:50:50,724 __bs: MainProcess Chudnovsky ... 132,000,000 iterations and 1550.18 seconds.
[DEBUG] 2022-10-03 09:51:12,742 __bs: MainProcess Chudnovsky ... 133,000,000 iterations and 1572.20 seconds.
[DEBUG] 2022-10-03 09:51:19,799 __bs: MainProcess Chudnovsky ... 134,000,000 iterations and 1579.26 seconds.
[DEBUG] 2022-10-03 09:51:28,824 __bs: MainProcess Chudnovsky ... 135,000,000 iterations and 1588.28 seconds.
[DEBUG] 2022-10-03 09:51:35,324 __bs: MainProcess Chudnovsky ... 136,000,000 iterations and 1594.78 seconds.
[DEBUG] 2022-10-03 09:51:48,419 __bs: MainProcess Chudnovsky ... 137,000,000 iterations and 1607.88 seconds.
[DEBUG] 2022-10-03 09:51:55,634 __bs: MainProcess Chudnovsky ... 138,000,000 iterations and 1615.09 seconds.
[DEBUG] 2022-10-03 09:52:04,435 __bs: MainProcess Chudnovsky ... 139,000,000 iterations and 1623.89 seconds.
[DEBUG] 2022-10-03 09:52:11,583 __bs: MainProcess Chudnovsky ... 140,000,000 iterations and 1631.04 seconds.
[DEBUG] 2022-10-03 09:52:17,222 __bs: MainProcess Chudnovsky ... 141,000,000 iterations and 1636.68 seconds.
[DEBUG] 2022-10-03 10:02:43,939 compute: MainProcess    Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 calulation Done! 141,027,339 iterations and 2263.39 seconds.
[INFO] 2022-10-03 10:09:07,119 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-03 10:09:07,119 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 10     Chudnovsky brothers  1988 
    π = (Q(0, N) / 12T(0, N) + 12AQ(0, N))**(C**(3/2))
 
[INFO] 2022-10-03 10:09:07,120 <module>: MainProcess Calculation took 141,027,339 iterations and 0:44:06.398345.

math_pi。Pi (b = 1000000) 快到一百万。大约快40倍。但它不能达到十亿,一百万是最多的数字。

GMPY内置看起来像:

python pi-pourri.py -v -d 1,000,000,000 -a 11
[INFO] 2022-10-03 14:33:34,729 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-03 14:33:34,729 compute: MainProcess Starting   const_pi() function from the gmpy2 library formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-03 15:46:46,575 compute: MainProcess    const_pi() function from the gmpy2 library calulation Done! 1 iterations and 4391.85 seconds.
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 11     const_pi() function from the gmpy2 library 
[INFO] 2022-10-03 15:46:46,575 <module>: MainProcess Calculation took 1 iterations and 1:13:11.845652.

内置的MPmath几乎一样快。慢约12%(6分钟):

python pi-pourri.py -v -a 12 -d 1,000,000,000  
[INFO] 2022-10-04 09:10:37,085 <module>: MainProcess Computing π to 1,000,000,000 digits.
[DEBUG] 2022-10-04 09:10:37,085 compute: MainProcess Starting   mp.pi() function from the mpmath library formula to 1,000,000,000 decimal places
[DEBUG] 2022-10-04 10:01:25,321 compute: MainProcess    mp.pi() function from the mpmath library calulation Done! 1 iterations and 3048.22 seconds.
[INFO] 2022-10-04 10:01:25,338 <module>: MainProcess Last 5 digits of π were 45519 as expected at offset 999,999,995
[INFO] 2022-10-04 10:01:25,340 <module>: MainProcess Calculated π to 1,000,000,000 digits using a formula of:
 12     mp.pi() function from the mpmath library 
[INFO] 2022-10-04 10:01:25,343 <module>: MainProcess Calculation took 1 iterations and 0:50:48.250337.

正好是3![弗林克教授(辛普森一家)]

开玩笑,但这里有一个在c#(。微软网络框架。

using System;
using System.Text;

class Program {
    static void Main(string[] args) {
        int Digits = 100;

        BigNumber x = new BigNumber(Digits);
        BigNumber y = new BigNumber(Digits);
        x.ArcTan(16, 5);
        y.ArcTan(4, 239);
        x.Subtract(y);
        string pi = x.ToString();
        Console.WriteLine(pi);
    }
}

public class BigNumber {
    private UInt32[] number;
    private int size;
    private int maxDigits;

    public BigNumber(int maxDigits) {
        this.maxDigits = maxDigits;
        this.size = (int)Math.Ceiling((float)maxDigits * 0.104) + 2;
        number = new UInt32[size];
    }
    public BigNumber(int maxDigits, UInt32 intPart)
        : this(maxDigits) {
        number[0] = intPart;
        for (int i = 1; i < size; i++) {
            number[i] = 0;
        }
    }
    private void VerifySameSize(BigNumber value) {
        if (Object.ReferenceEquals(this, value))
            throw new Exception("BigNumbers cannot operate on themselves");
        if (value.size != this.size)
            throw new Exception("BigNumbers must have the same size");
    }

    public void Add(BigNumber value) {
        VerifySameSize(value);

        int index = size - 1;
        while (index >= 0 && value.number[index] == 0)
            index--;

        UInt32 carry = 0;
        while (index >= 0) {
            UInt64 result = (UInt64)number[index] +
                            value.number[index] + carry;
            number[index] = (UInt32)result;
            if (result >= 0x100000000U)
                carry = 1;
            else
                carry = 0;
            index--;
        }
    }
    public void Subtract(BigNumber value) {
        VerifySameSize(value);

        int index = size - 1;
        while (index >= 0 && value.number[index] == 0)
            index--;

        UInt32 borrow = 0;
        while (index >= 0) {
            UInt64 result = 0x100000000U + (UInt64)number[index] -
                            value.number[index] - borrow;
            number[index] = (UInt32)result;
            if (result >= 0x100000000U)
                borrow = 0;
            else
                borrow = 1;
            index--;
        }
    }
    public void Multiply(UInt32 value) {
        int index = size - 1;
        while (index >= 0 && number[index] == 0)
            index--;

        UInt32 carry = 0;
        while (index >= 0) {
            UInt64 result = (UInt64)number[index] * value + carry;
            number[index] = (UInt32)result;
            carry = (UInt32)(result >> 32);
            index--;
        }
    }
    public void Divide(UInt32 value) {
        int index = 0;
        while (index < size && number[index] == 0)
            index++;

        UInt32 carry = 0;
        while (index < size) {
            UInt64 result = number[index] + ((UInt64)carry << 32);
            number[index] = (UInt32)(result / (UInt64)value);
            carry = (UInt32)(result % (UInt64)value);
            index++;
        }
    }
    public void Assign(BigNumber value) {
        VerifySameSize(value);
        for (int i = 0; i < size; i++) {
            number[i] = value.number[i];
        }
    }

    public override string ToString() {
        BigNumber temp = new BigNumber(maxDigits);
        temp.Assign(this);

        StringBuilder sb = new StringBuilder();
        sb.Append(temp.number[0]);
        sb.Append(System.Globalization.CultureInfo.CurrentCulture.NumberFormat.CurrencyDecimalSeparator);

        int digitCount = 0;
        while (digitCount < maxDigits) {
            temp.number[0] = 0;
            temp.Multiply(100000);
            sb.AppendFormat("{0:D5}", temp.number[0]);
            digitCount += 5;
        }

        return sb.ToString();
    }
    public bool IsZero() {
        foreach (UInt32 item in number) {
            if (item != 0)
                return false;
        }
        return true;
    }

    public void ArcTan(UInt32 multiplicand, UInt32 reciprocal) {
        BigNumber X = new BigNumber(maxDigits, multiplicand);
        X.Divide(reciprocal);
        reciprocal *= reciprocal;

        this.Assign(X);

        BigNumber term = new BigNumber(maxDigits);
        UInt32 divisor = 1;
        bool subtractTerm = true;
        while (true) {
            X.Divide(reciprocal);
            term.Assign(X);
            divisor += 2;
            term.Divide(divisor);
            if (term.IsZero())
                break;

            if (subtractTerm)
                this.Subtract(term);
            else
                this.Add(term);
            subtractTerm = !subtractTerm;
        }
    }
}