自动搬运

查看原文

来自洛谷，原作者为

	TBB_Nozomi **

搬运于2025-08-24 21:28:24，当前版本为作者最后更新于2020-06-14 01:58:14，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

---

前记：这道题可能是全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})$，有：

$$\boxed{|R_n(x)|=\frac{\mathrm{e}^{\theta x}x^{n+1}}{(n+1)!}\leq \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}}$。