题解 P1729 【计算e】

· · 题解

前记：这道题可能是全Luogu唯一的一道又有高精度又有\mathrm{e}的题目，我就用这道题来练手了(

这道题的普遍做法（也是较为常见的正解）是通过提前打表存入程序里，再根据输入的参数做格式化输出。这题也不是什么考场上的题，因此这么做其实相当的靠谱。一种常见的方式就是在Wolfram Mathematica里面输入N[\mathrm{e},10001]，然后你就可以立即得到一份小数位10000位的\mathrm{e}，代入程序即可，这里不再赘述。这里讲述两种非打表的高精度解法。

由题目用意，自然对数的底数\mathrm{e}的其中一种定义为:

\mathrm{e}=\sum\limits_{n=0}^{\infty}\frac{1}{n!}

那么我们的两种做法就围绕着这个式子展开。

做法一

通分+高精度除法

前置技能：简单的微积分（泰勒展开或各种中值定理）或初高中阶段简单的不等式知识；P5432 高精度除法。

观察定义给出的这个式子，我们容易发现n不一定非得算到无穷大，我们只需要使用前面的若干位就可以得到\mathrm{e}的近似值。那么我们要取多少位呢？高数/数分里面，关于Taylor展开的前若干项的误差，有一个Lagrange余项的估计：

\mathrm{e}_ {N}= \sum\limits_{k=0}^{N}\frac{1}{k!}\qquad R_N=\mathrm{e}-\mathrm e_N=\frac{e^\theta}{(N+1)!}

其中\theta\in(0,1)。即便你并不知道什么是拉格朗日余项，也可以通过不等式的方式分析误差：考虑余项中的每一项\frac{1}{(k+1)!}的分母，我们只取前N个数和最后的两个数，那么分式将小于\frac{1}{N!k(k+1)}=\frac{1}{N!}(\frac{1}{k}-\frac{1}{k+1})，因此余项有：

R_N=\sum\limits_{k=N+1}^{\infty}\frac{1}{k!}<\frac{1}{N!}(\frac{1}{N+1}+\sum\limits_{k=N+1}^{\infty}(\frac{1}{k}-\frac{1}{k+1}))<\frac{1}{N!}

假定本题要求输出的相对误差\delta<10^{-k}（本题的输出结果是一个固定的值，很容易把绝对误差变成相对误差），那么有

\boxed{\delta_N=\frac{R_N}{\mathrm{e}}<\frac{1}{N!}<10^{-k}=\delta}

在本题中有k\leq10000，因此可以得到N\leq3249。考虑到一些误差因素，可以适当的将N放大一下以得到精度更高的结果。

那么，就按照这一结果计算就完事了吗？鉴于高精度小数的常数，可能对比较大的数据会有些困难。考虑到本题中需要精度位数为k，照公式计算中的每一位都应该至少也要有k位。由Stirling公式我们可以大概的估计N与k有关系k\sim N\log N，而逐项使用高精度定点数与单精度整数的除法的时间复杂度为\mathcal{O}(k)（高精度浮点数的常数只会更大），则总时间复杂度为\mathcal{O}(kN)=\mathcal{O}(\frac{k^2}{\log N})，对于k\sim 10000部分而言可能会比较危险，有TLE的可能。

怎么样优化掉这个\mathcal{O}(nk)的累加和？其实前面的Qingyu大佬有讲过，我们可以把所有的这些分式通分，变成：

\mathrm{e}_ {N}= \sum\limits_{k=0}^{N}\frac{1}{k!}=\frac{1}{N!}\sum\limits_{k=0}^{N}A_N^k=\frac{1}{N!}\sum\limits_{k=0}^{N}\prod\limits_{t=N-k+1}^{N}t

这样，我们就可以分别计算分子和分母，再通过高精度除法得到最后的结果，这一过程中仅有最后的高精度除法是稍有误差的。前面使用高精度整数加法和与单精度的乘法，其时间复杂度大致与前面的以一致，为较小常数的\mathcal{O}(\frac{k^2}{\log N})；后面如果采取高精度除法而非模拟试除法，则可以优化到\mathcal{O}(k\log k)。经实验该方法本题可过。

本方法部分代码：

int main()  {
    int k;
    cin>>k;
    tbb::_LFloat_prec= (k/4)+2; //万进位高精
    LInt S= 1, P= 1, T= 1;
    int N;
    for(N=1; T.digit()<=tbb::_LFloat_prec*4; )  T*= (++N); //先估计分母阶乘的位数
    for(int i=N; i>0; i--)  {
        P*= i;  S+= P;
    }
    LFloat F=S; F/=T; //做高精度浮点数除法
    string ans= F.print_str();
    const char* out= ans.c_str();
    putchar('2');   putchar('.');   putchar('\n');
    for(int T=1; T<=k; ++T) {
        putchar(out[1+T]);
        if(T%50==0) putchar('\n');
        else    if(T%10==0) putchar(' ');
    }
    return 0;
}

做法二

实现高精度浮点数的exp函数

前置技能：比较熟练的微积分；《理性愉悦：高精度数值计算》；高精度浮点数的加减乘法，高精度小数的正整数幂。（不需要知道高精度除法真高兴）

既然本题是Luogu罕见的高精度又带\mathrm{e}，那自然就不会放过这个机会测试函数\exp啦。我们可以考虑实现\exp(x)=\mathrm{e}^x，再进行计算\exp(1)即可。那么我们接下来将着手分析并实现\exp。

根据泰勒展开，我们可以将函数\mathrm{e}^x展开成幂级数形式，并使用拉格朗日余项，我们有：

\mathrm{e}^x=\sum\limits_{k=0}^{\infty}\frac{x^k}{k!}= \sum\limits_{k=0}^{N}\frac{x^k}{k!}+R_N(x)\quad |R_N(x)|=\frac{\mathrm{e}^{\theta x}x^{N+1}}{(N+1)!}\leq\frac{\mathrm{e}^x}{(N+1)!}

其中\theta\in(0,1)。由此得知如此展开以后的相对误差也不会超过\frac{1}{N!}。

同上面一样，虽然看起来本题可以计算了，然而仔细分析一下，高精度浮点数的乘法次数至少是\mathcal{O}(N)的，而高精度乘法的时间复杂度是\mathcal{O}(k\log k)的，放一起估计一下也有\mathcal{O}(k^2\log k)，那肯定会超时。虽然针对本题，后面会讲到如何优化掉这部分乘法，但总的来说还是要尽可能的降低这部分乘法的次数。

虽然与本题无关，但是这里还是简单的提一句：对于整个实数范围内的\exp求值，我们可以将它规约到[0,1]的范围内。对于x<0的，有\mathrm{e}^{x}=1/\mathrm{e}^{-x}可以化为x>0的求值；对于x>1的，由常见的浮点数存储格式，可以认为x=x_0\times 10^E（其中x_0\in (0,1)并使用高精度整数来存储）(P.S. 这很像但并不是科学记数法，科学计数法可以规约到这种形式)，那么有\exp(x)=\exp(x_0)^{10^E}，可以化为(0,1)范围内的求值再做个高精度浮点数的正整数次幂得到。

参照PPT，我们可以进一步的将x的范围从(0,1]放缩到(0,10^{\lfloor-\sqrt{k}\rfloor})，然后仿照上面讨论x>1的部分那样提出10^{\sqrt{x}}放到\exp的外面再做快速幂。对于x\in(0,10^{\lfloor-\sqrt{k}\rfloor})，有：

\mathrm{e}^x10^{-(n+1)\lfloor\sqrt{k}\rfloor}\quad\delta=\frac{R_n(x)}{\mathrm{e}^x}<10^{-(n+1)\lfloor\sqrt{k}\rfloor}}

则只需要取n=\mathcal{O}(\sqrt{k})，无论是\exp的里面还是外面都只需要\mathcal{O}(\sqrt{k})次长度为k的高精度乘法，因此本方法的总的时间复杂度为\mathcal{O}(k^{1.5}\log k)。

虽然时间复杂度理论上满足了，然而实践操作中，如果是万进压位的高精度浮点数，对于k=10000的点，需要做大概500次的长度为4096~8192位的高精度乘法，这部分对常数的要求很大。我们可以采取一些比较通用的办法来进行优化。比如比较的经典的包括预处理循环卷积的单位根；又或者是对为循环卷积弄个缓存并适当调整乘法的位置，以保证不总是重复对同一个数组做DFT。当然也可以调整n=\alpha\sqrt{k}中的常数\alpha，平衡放缩前与放缩后的乘法次数。经过这些优化以后，我针对本题可以做到平均1.2s一个点，也许可以有更快的常数。

然而，针对本题，有一种更好的优化方法。如果你采用的是形如x=x_0\times 10^E这样的标准的存储法，你会发现无论是放缩前还是放缩后，x中的有效数位其实都只有一位。原因是不言自明的。因此在累加累乘的时候，只要在进行乘法时先判断两边的有效数位，其中的一边大概小于另一边的对数的级别的时候，就可以直接用朴素的按位数乘法累乘而不用经过循环卷积。这样可以减少至少一半的循环卷积次数，大大的提高运行效率。

部分核心代码如下：

LFloat pow_pow10(const LFloat & x, int m)   {
    LFloat S= x, S_2, S_3;
    for(int i=0; i<m; ++i)  {
        S= S.pow2();
        S_2= S*S;   S_3= S_2* S; S= S_2* S_3;
    }
    return S;
}// return x^(10^m)
LFloat exp(const LFloat& x)   {
    if(x.isNaN())   return x;
    if(x.isinf() && x.positive())   return x;
    if(x.isinf() && x.negative())   return 0;
    if(x==0)    return 1;
    if(x<0)     return 1/exp(-x);
    if(x > (double(_LFloat_prec)+INT_MAX)*std::log(10000))   return LInt("inf");
    if(x < (double(INT_MIN))*std::log(10000))    return 0;

    int precision= _LFloat_prec * 4;
    if(x>1) {
        _LFloat_prec= (precision + 4 * (x.pow+x.base.d) )/4 + 1;
         LFloat res= pow_pow10( exp(LFloat(x.base, -x.base.d)), 4 * (x.pow+x.base.d) );
         _LFloat_prec= precision / 4;
         return res;
     }
     int bound= int(std::sqrt(precision)/2.3);
     LFloat B= pow10<LFloat>(-bound);
     if(x>B) {
          int delta_pow= 4*x.pow + x.base.digit() + bound;
          _LFloat_prec= (precision + delta_pow)/4 + 1;
          LFloat x_= mul_pow10(x, -delta_pow);
          LFloat res= pow_pow10(exp(x_), delta_pow);
          _LFloat_prec= precision / 4;
          return res;
      }

      _LFloat_prec= (precision + Log_2(precision))/4 + 1;
      int N= precision / bound +1;
      LFloat S= 1;
      for(int i=N; i>0; --i)  S= S/i*x + 1;
      _LFloat_prec= precision/4;
      S.sho();
      return S;
}

int main()  {
    int k;
    cin>>k;
    tbb::_LFloat_prec= (k/4)+2;
    LFloat F= exp(LFloat(1));
    string ans= F.print_str();
    const char* out= ans.c_str();
    putchar('2');   putchar('.');   putchar('\n');
    for(int T=1; T<=k; ++T) {
        putchar(out[1+T]);
        if(T%50==0) putchar('\n');
        else    if(T%10==0) putchar(' ');
    }
    return 0;
}

其它代码注释

LInt表示高精度整数；LFloat表示高精度浮点数，由一个LInt的base和一个int的pow组成，其表示的数为\mathrm{base}\times10000^{\mathrm{pow}}。