高精度计算pi [P1727]

TBB_Nozomi · 2021-04-30 04:04:07 · 题解

blog:Masof TBB

前记：本题是目前洛谷所有含有“高精度”的 tag 中唯一一个含有高精度圆周率计算的题目。如果对本题做出扩展，也可以将本题当作高精度反三角函数的题目来做。因此我就来这个题下手了（

在本文中，本文出示本题解答的四种方法，并简要证明其它题解中出现但尚未证明的一些恒等式，帮助各位读者理解。在不另外加以说明的情况下，本题默认 N 表示题目所需的有效数位（即精度为 10^{-N}），n 表示迭代、累加所需要的项的个数。

0. 打表

前人之述备矣。如果你想用打表来求解此题，唯一需要注意的问题是将内容格式化为程序所能接受的形式。\pi 的数据的来源可以是网上搜索，也可以是线下自己手写高精度计算（如果能读者做到这一点，想必是不需要本篇题解也能通过此题的了）。一种来源是在软件 Mathematica 上输入 N[Pi, 10001] 即可得到满足题目需求的圆周率小数点后 10000 位的全文本。

1. 级数与迭代

前置知识：一个你觉得好用的好计算的 \pi 的求和式或者迭代式。可选的有简单的微积分知识（泰勒展开和各种中值定理）用于得到计算式和估计尾项；各种高中的三角函数公式；高精度浮点数除法，或者是一点简单的高精度定点数对低精度数的乘法和除法。

和 \mathrm{e} 不一样，\pi 的计算由来已久，在数学史上发展的几乎每一个阶段都会生成出各种各样求 \pi 的公式。如果不追求效率和精度，你甚至可以使用蒙特卡洛法来估算 \pi。下面给出一个不适宜本题使用的级数的例子：一个最有名的 \pi 的级数表达式是

称为 Leibniz-Gregory 级数。这一式子本质上来源于反正切函数 \arctan x 在 x=0 处的 Maclaurin 展开式：

然后将 x=1 的值代入。虽然这个式子相当经典，但是这个级数的收敛速度很慢，如果做一个估计的话，会发现它的余项大概在 \mathcal{O}(\frac{1}{n}) 这个量级（其中 n 为求和所使用的项数）。具体的误差分析放到后面高精度反正切函数的分析部分，但可以确定它并不适宜来在限定时间中做出此题。

从 \arctan 1 出发，利用正切函数相关的公式，如 \tan (u+v) = \frac{\tan{u}+\tan{v}}{1-\tan{u}\tan{v}}，我们可以将一个反正切函数，写成多个反正切函数的和与差。举个例子，如果我们令 u=\arctan\frac{1}{2}, v=\arctan\frac{1}{3}，代入上面的和角公式，就可以得到 u+v=\arctan 1。同理，取 u=\arctan\frac{1}{5} 和 v=\arctan\frac{1}{239}，则可以得到 4u-v=\arctan 1。反复使用公式，你甚至能得到 \pi=32\arctan\frac{1}{10}-4\arctan\frac{1}{239}-16\arctan\frac{1}{515} 这样看起来非常壮观的公式。这类 Machin 公式的一个共同思想，就是利用正切函数的角加减变换，将 \arctan 1 变成若干个 \arctan x（其中 x<1）的和，使得反正切的麦克劳林级数收敛速度更快。为方便计算，通常都会把分解后的 x 写成为分子为 1 的单位分数。

我们取其中一个公式作为例子。由 \arctan\frac{1}{2}+\arctan\frac{1}{3}=\arctan 1，我们有：

哪怕仅仅只是从数列的项的大小来看，我们也能预感到这个级数的收敛速度应该远快于最开始的那个例子。事实上，取 \arctan\frac{1}{k} 带入到麦克劳林展开式里面，则它的余项大概在 \mathcal{O}(\frac{1}{k^N}) 这一量级。这样的公式已经可以用来做本题了。

在一些人题解中，提到了下面这个公式

其中 (2k+1)!! = 1\times3\times5\times\dots\times(2k+1) 表示双阶乘。这也是一个收敛速度比较快的级数公式。根据 Pi Formula 页面上的信息，这一公式每计算一项就能提供约 2 bits 的信息，也就是平均计算 3 项可以使你十进制的精确度增加一位。相关的证明以及误差分析将放在第2.5节说明。

上面的公式都有一个共同点：他们的数列的通项形式上都相对的简单，而且普遍的能使用前一项去推出后一项。利用这一特点，我们可以调整项的计算方式，来得到更加简洁，且精确度更高计算过程。以上面的 (Euler-T) 的式子为例，假设我们已经确定使用前 N 项求和（k=0 的项视为第 0 项），则可以调整计算过程如下：

由这个式子，我们从 S=1 出发，倒着枚举 k 从 N 到 1，每一步都取 1+\frac{k}{2k+1}S \to S 的赋值操作。这样只需要我们有一个高精度对低精度乘除的操作就行。不然你可能需要通分，然后再做一个高精度除法。

时间复杂度：设所需的精度为 N 位，上述方法通常的给出了一个需要做 \mathcal{O}(N) 次迭代操作（这是因为公式中每多计算一项，都只能提供常数的有效数位），每一次迭代操作都需要 \mathcal{O}(N) 的计算量（用于加法、对低精度的乘除法）。综合起来总的时间复杂度为 \mathcal{O}(N^2)，在本题的数据范围内应该可以全部通过。

部分代码样例：（使用 (4) 式子，且采取先通分再高精度除法的方法）

LFloat _pi(void)    {
    static LFloat inner_pi= 3;
    static int pi_precsion= 0;
    if(_LFloat_prec < pi_precsion)  { //缓存pi的计算结果，若精度满足要求则直接输出
        LFloat ret = inner_pi;
        ret.sho();  return ret;
    }
    pi_precsion= _LFloat_prec;
    if(_LFloat_prec > 1000) return inner_pi = 4 * arctan(1);
    static LInt numerator= 1, denominator= 1; // 缓存上次计算时的分子和分母
    static LInt iter_factorial= 1;
    static int done_iter_times= 0;
    int new_iter_times= 6.09832411290916 + 13.28771237954945 *_LFloat_prec;
    LInt delta_numerator=0, delta_denominator= 1, delta_iter_factorial= 1;
    for(int k= new_iter_times; k > done_iter_times; --k)    {   // 计算新增加的项的部分的分子和分母
        delta_numerator+=   delta_denominator;
        delta_numerator*=   k;
        delta_denominator*= 2*k+1;
        delta_iter_factorial*=  k;
    }
    numerator= numerator * delta_denominator + iter_factorial * delta_numerator;
    denominator*= delta_denominator; // 计算通分后的新的分子和分母
    iter_factorial*= delta_iter_factorial;
    return Div_LInt(numerator, denominator) * 2;
}

2. 高精度反正切函数

前置知识：基本的微积分知识，时间复杂度为 \mathcal{O}(N \log N) 的高精度乘法(P1919)。可选的其它知识与材料：2012 年 WC 课件《理性愉悦——高精度数值计算》；时间复杂度为 \mathcal{O}(N \log N) 的高精度浮点数除法(P5432)和开平方根(P2293)。

要实现高精度的反正切函数，最简单的方法先将 x 规约到一个适合计算的范围内，然后在 x=0 处做 Taylor 展开，利用 Maclaurin 公式，取前面的若干项做求和，并估计误差项使其满足精度要求。由函数的奇偶性，x<0 的部分是可以转化为 x>0 的。

首先我们需要对 \arctan x 高次求导，由归纳法或者是查表有

由此可以得到 \arctan^{(2k+1)}0 = (-1)^k(2k)! 和 \arctan^{(2k)}0 = 0，代入 Taylor 展开式即有 (2) 式。

我们将相邻的正负符号相反的两项看作一组项，现在我们仅仅只取 (2) 式中幂指数小于等于 4n-1 的项，也即选取里面的前 2n 个项求和，用来作为对 \arctan x 的估计，则余项可以被写作为 R_{4n-1}(x)。根据 Taylor 展开的 Lagrange 余项，有

如果我们选取 x>0，则有放缩估计：

这意味着我们能估计的相对精度 \delta_{4n-1}=\frac{R_{4n-1}(x)}{\arctan x} 的上界将取决于我们选取的 x 和 n 的大小。如果我们的 x 的界仅仅能达到 1，则这种方法是不能拿来估计我们所需要的计算的项的。事实上，当 x=1 时，尾项有

相比于这题所需要的精度而言，(1) 的收敛速度太慢了。这启发我们（如同课件所说的那样）先将用于展开计算 x 的数量级变得比较小，再进行计算和转化。

利用正切函数的倍角公式，我们可以得到 \arctan x= 2\arctan f(x)，其中

其中当 x>1 时 f(x)<1。而 x\leq 1 时 f(x)\leq \frac{x}{2}，因此可以认为它跟 exp 函数一样，拥有一个等比例的放缩函数，且每放缩一次所需的计算量为 \mathcal{O}(N \log N)（计算量源自于放缩函数中的开平方根和除法运算）。

另一方面，放缩以后利用相邻项之间的递推，每计算一项就会进行一次高精度乘法。不妨设最后我们将 x 放缩到 [0,10^{-B}) 的范围内，则 Taylor 展开计算的项数 n 应该满足 10^{-4nB}\leq10^{-N} 即 n\geq\frac{N}{4B}。在放缩上，将 [0,1) 放缩到 [0,10^{-B}) 上需要的放缩次数 k 应该满足 2^{-k}\leq 10^{-B} 即 k\geq \log_2{10}B。则本函数的总的时间复杂度为 \mathcal{O}((n+k)N\log N)，由基本不等式，应该选取 B=\mathcal{O}(\sqrt{N})，可以使算法的复杂度降低为 \mathcal{O}(N^{1.5}\log N)。注意到放缩的函数中，无论是开根号还是除法，其时间复杂度虽然与乘法相同，但常数却几倍于乘法，因此应该适当减小 B 的值。

部分代码样例如下。经测试，在 -O2 情况下可以至多 2.96s 的时间通过此题。

LFloat arctan(const LFloat& x)  {
    if(x.isNaN()||x.zero())   return x;
    if(x.isinf() && x.negative())   return arctan(-1)*2;
    if(x.isinf() && x.positive())   return arctan(1)*2;
    if(x<0) return -arctan(-x);
    if(x>1) return arctan((sqrt(x*x+1)-1)/x)*2; // 规约x到[0,1]中

    struct{
        LFloat operator()(const LFloat& t) {return t / (sqrt(t*t+1)+1);}
    } scale_func; // 放缩函数
    int precision= _LFloat_prec * 4;
    int bound= int(std::sqrt(0.2006866637759875 * precision / 16));    //bound= Sqrt(2 Lg 2*p/3)/16
    LFloat B= pow10<LFloat>(-bound);
    int n= precision/bound +1;  //expansion terms count
    int k= 0;   //scale times

    _LFloat_prec= (precision + Log_2(precision) + 3.322*bound)/4 + 1;
    LFloat x_scaled= x;
    x_scaled.sho();
    while(x_scaled > B) {   //scaling
        x_scaled= scale_func(x_scaled);
        ++k;
    }
    LFloat y_scaled= 0, x2= x_scaled * x_scaled;
    for(int i= 4*n-1; i>=1; i-=2)   y_scaled= -x2*y_scaled + LFloat(1.0)/i; // 求较小的x的近似值
    y_scaled*= x_scaled;
    LFloat& y= y_scaled;
    for(int i= 0; i<k; ++i) y= y*2;
    _LFloat_prec= precision/4;
    y.sho();
    return y;
}

2.5. 欧拉变换与超几何函数

阅读前提醒：本部分用于加速反正切函数的计算，但是不能改善其时间复杂度。同时，本部分也会介绍广受好评的公式 (4) 的来源。跳过本部分不会影响本题的解答。

前置知识：第2节前置知识；组合数相关知识或生成函数。可选其它补充知识：超几何函数（《具体数学》）。相关的资料可以参考 Wolfram Mathworld 的 Euler Transform。

在继续前进之前，我们先来了解一个恒等式。假设 S=a_1-a_2+a_3-a_4+a_5-\dots 是一个收敛的交错级数，定义

则 (s_n)\triangleq(s_{n,1})_{n=1}^{\infty} 就构成了一个数列。则我们有

后面的级数将与前一个级数收敛到同一个数，且收敛速度更快。这被称作欧拉改进收敛法，也被称为欧拉变换。这里用生成函数简单证明一下这个式子的收敛性：在定义 (11) 中，利用组合数，可知 s_n 中含有的 a_k 的系数为 (-1)^{k-1}\binom{n-1}{k-1}。因此改进后的级数中，所含有的 a_k 的系数为

其中 [x^k]f(x) 表示 f 的展开式中 x^k 项的系数。上式表明改进后的收敛级数在形式上与原先的级数具有相同的组成部分，因此它们（如果收敛）收敛于同一个值。

接下来我们看看这个恒等式能用在什么地方。仔细一看，(1) 刚好是一个交错级数。取 a_k=\frac{1}{2k-1}，由归纳法可以计算出来 s_n=\frac{2^{n-1}(n-1)!}{(2n-1)!!}，再代回恒等式的右端，它正好就是式子 (4) 两端同时除以 2 的结果！这样，我们就严格的推导出了 (4)。

我们再来看一下这个恒等式能有什么用途。观察 (2) 发现它也是一个交错级数，但是直接运用到恒等式上面未免过于繁琐。这里，我们引入超几何函数 {}_2\!F_1，它被定义为

其中，打括号的幂指数部分表示上升幂(不知道为什么洛谷突然不支持 \overline 指令了) 。特别的，有

不加证明的，这里给出一个超几何变换（一个恒等式）:

这个变换公式被称为反射定律，也被称为欧拉超几何变换。当 (15) 中取 z=-1 时就构成了欧拉改进收敛法的一个特例。现在将这个公式直接套用在 (14) 的左边，就可以得到一个新的 \arctan x 的表达式：

这个式子已经变成了一个计算量常数相对小的正项级数了，可以跑出一个相当不错的时间。按照第2节内的方法调整放缩和展开系数，可以达到常数比较小的 \mathcal{O}(N^{1.5}\log N) 的时间复杂度。

部分代码如下。在开启 -O2 选项的时候成功跑进了 1s 内，不开启时最大的点也只有 4.90s。

LFloat arctan2(const LFloat& x)  {
    if(x.isNaN()||x.zero())   return x;
    if(x.isinf() && x.negative())   return arctan(-1)*2;
    if(x.isinf() && x.positive())   return arctan(1)*2;
    if(x<0) return -arctan(-x);
    if(x>1) return arctan((sqrt(x*x+1)-1)/x)*2;

    struct{
        LFloat operator()(const LFloat& t) {return t / (sqrt(t*t+1)+1);}
    } scale_func;
    int precision= _LFloat_prec * 4;
    int n= std::ceil(std::sqrt(precision/2)*12) + 1;
    int bound= std::ceil(1.0*precision/2/n) + 1;  //expansion terms count
    LFloat B= pow10<LFloat>(-bound);
    int k= 0;   //scale times

    _LFloat_prec= (precision + Log_2(precision) + 3.322*bound)/4 + 1;
    LFloat x_scaled= x;
    x_scaled.sho();
    while(x_scaled > B) {   //scaling
        x_scaled= scale_func(x_scaled);
        ++k;
    }
    LFloat z= 1 / (1 + x_scaled * x_scaled);
    LFloat y_scaled= 1, t= 1-z;
    for(int i= n; i>= 1; --i)   y_scaled= y_scaled*t*2*i/(2*i+1) + 1;
    y_scaled*= x_scaled * z;
    LFloat& y= y_scaled;
    for(int i= 0; i<k; ++i) y= y*2;
    _LFloat_prec= precision/4;
    y.sho();
    return y;
}

3. AGM 方法

前置知识：《理性愉悦——高精度数值计算》。

参考论文：Eugene Salamin 于 1976 年发表于 Mathematics of Computation 上的论文 Computation of pi Using Arithmetic-Geometric Mean。

本来这篇题解到这里就该结束了，但没想到课件的后面提出另一个作为目前计算初等函数时间复杂度最低的方法：AGM 方法。囿于这部分相关的数学知识不熟悉，本题解无法提供该方法正确性的证明。

下面简单的介绍一下这个方法。假设给出两个正数 a 和 b，取两个序列 (a_n) 和 (b_n)，其中 a_0=a, b_0=b，按下列方法做递推：

这两个序列的取值分别为上一项的两个数的算数平均值和几何平均值序列。容易证明这两个序列一个单调递增一个单调递减，且对应项的差会趋近于 0。因此这两个序列会趋近于同一个值，称为两个数的算数几何平均值，记为 \mathrm{agm}(a,b)，这个值等于一个椭圆积分的表达式。1975 年，Eugene Salamin 在他的论文里提出，如果取得适当的 a 和 b，则它们的 \mathrm{agm} 值会变成一个特殊的数，并提出当 a=1 且 b= \frac{1}{\sqrt{2}} 时，有

这个 \pi 的逼近数列的误差满足

可以发现 \pi_n 的精确度是指数级别增长的。取 n=\left\lceil\log_2N\right\rceil，当 N>1 时有 -\lg\left|\pi-\pi_n\right|>N，此时迭代即可达到精度要求。

时间复杂度：每迭代一项需要做一次高精度加法、一次高精度乘法和一次高精度开平方根，也就是 \mathcal{O}(N \log N)。n=\mathcal{O}(\log N)，因此总的时间复杂度为 \mathcal{O}(N \log^2N)，确实很快。

部分代码如下。在不开启 -O2 的情况下，最慢的点需要 1.2s，若开启则可以加速到 300ms 以内。

LFloat agm_pi(void) {
    int precision= tbb::LFloat::precision() * 4;
    tbb::LFloat::precision(std::ceil((precision+Log_2(precision))/4.0)); // 调整过程精度
    LFloat a= 1, b= 1/sqrt(LFloat(2));
    int n= Log_2(precision);
    LFloat S= pow(b-a, 2);
    for(int i=1; i<=n; ++i) {
        LFloat an= (a+b)/2, bn= sqrt(a*b);
        S+= (1<<i)*pow(bn-an, 2);
        a= an, b= bn;
    }
    LFloat ans= (4*a*a)/(1-S);
    tbb::LFloat::precision(precision/4);
    ans.sho();
    return ans;
}