计算e
题目背景
《爱与愁的故事第二弹·compute》最终章。
自然对数的底数 e e e 是一个著名的无理数，其近似值为 2.718281828 ⋯ 2.718281828\cdots 2.718281828⋯ 有计算 e e e 的公式如下：
e = ∑ n = 0 ∞ 1 n ! e=\sum_{n=0}^{\infty}\frac{1}{n!} e=n=0∑∞n!1
其中 n ! n! n! 表示 n n n 的阶乘，即 n ! = 1 × 2 × 3 × ⋯ × n n!=1\times 2\times 3\times \cdots \times n n!=1×2×3×⋯×n。
题目描述
月落乌啼竟然这么快就回复了圆周率小数点后10000位？！不可能，他肯定求了别人。爱与愁大神再次为难月落乌啼：“帮我算一算 e e e 后 n ( n ≤ 10000 ) n(n \le 10000) n(n≤10000) 位，速度！！！”月落乌啼想求别人，结果他发现由于刚才跟你通话已经用完了手机的所有电。关键时刻只能靠自己。如果现在你是他，你会怎么编这个程序？
输入格式
只有一行： n n n
输出格式
若干行。
第一行： 2 2 2.
第二行开始： e e e 的小数部分。 10 10 10 个数一个空格， 50 50 50 个数一次回车。
样例 #1
样例输入 #1
100
样例输出 #1
2.
7182818284 5904523536 0287471352 6624977572 4709369995
9574966967 6277240766 3035354759 4571382178 5251664274
提示
30
%
30\%
30% 数据：
n
≤
1000
n \le 1000
n≤1000
100
%
100\%
100% 数据：
n
≤
10000
n \le 10000
n≤10000
时限：全部1秒
前记：这道题可能是全Luogu唯一的一道又有高精度又有e\mathrm{e}e的题目，我就用这道题来练手了(
这道题的普遍做法（也是较为常见的正解）是通过提前打表存入程序里，再根据输入的参数做格式化输出。这题也不是什么考场上的题，因此这么做其实相当的靠谱。一种常见的方式就是在Wolfram Mathematica里面输入N[e\mathrm{e}e,10001]，然后你就可以立即得到一份小数位10000位的e\mathrm{e}e，代入程序即可，这里不再赘述。这里讲述两种非打表的高精度解法。
由题目用意，自然对数的底数e\mathrm{e}e的其中一种定义为:
e=∑n=0∞1n!\mathrm{e}=\sum\limits_{n=0}^{\infty}\frac{1}{n!} e=n=0∑∞n!1
那么我们的两种做法就围绕着这个式子展开。
做法一
通分+高精度除法
前置技能：简单的微积分（泰勒展开或各种中值定理）或初高中阶段简单的不等式知识；P5432 高精度除法。
观察定义给出的这个式子，我们容易发现nnn不一定非得算到无穷大，我们只需要使用前面的若干位就可以得到e\mathrm{e}e的近似值。那么我们要取多少位呢？高数/数分里面，关于Taylor展开的前若干项的误差，有一个Lagrange余项的估计：
eN=∑k=0N1k!RN=e−eN=eθ(N+1)!\mathrm{e}_ {N}= \sum\limits_{k=0}^{N}\frac{1}{k!}\qquad R_N=\mathrm{e}-\mathrm e_N=\frac{e^\theta}{(N+1)!} eN=k=0∑Nk!1RN=e−eN=(N+1)!eθ
其中θ∈(0,1)\theta\in(0,1)θ∈(0,1)。即便你并不知道什么是拉格朗日余项，也可以通过不等式的方式分析误差：考虑余项中的每一项1(k+1)!\frac{1}{(k+1)!}(k+1)!1的分母，我们只取前N个数和最后的两个数，那么分式将小于1N!k(k+1)=1N!(1k−1k+1)\frac{1}{N!k(k+1)}=\frac{1}{N!}(\frac{1}{k}-\frac{1}{k+1})N!k(k+1)1=N!1(k1−k+11)，因此余项有：
RN=∑k=N+1∞1k!<1N!(1N+1+∑k=N+1∞(1k−1k+1))<1N!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!} RN=k=N+1∑∞k!1<N!1(N+11+k=N+1∑∞(k1−k+11))<N!1
假定本题要求输出的相对误差δ<10−k\delta<10^{-k}δ<10−k（本题的输出结果是一个固定的值，很容易把绝对误差变成相对误差），那么有
δN=RNe<1N!<10−k=δ\boxed{\delta_N=\frac{R_N}{\mathrm{e}}<\frac{1}{N!}<10^{-k}=\delta} δN=eRN<N!1<10−k=δ
在本题中有k≤10000k\leq10000k≤10000，因此可以得到N≤3249N\leq3249N≤3249。考虑到一些误差因素，可以适当的将NNN放大一下以得到精度更高的结果。
那么，就按照这一结果计算就完事了吗？鉴于高精度小数的常数，可能对比较大的数据会有些困难。考虑到本题中需要精度位数为kkk，照公式计算中的每一位都应该至少也要有kkk位。由Stirling公式我们可以大概的估计NNN与kkk有关系k∼NlogNk\sim N\log Nk∼NlogN，而逐项使用高精度定点数与单精度整数的除法的时间复杂度为O(k)\mathcal{O}(k)O(k)（高精度浮点数的常数只会更大），则总时间复杂度为O(kN)=O(k2logN)\mathcal{O}(kN)=\mathcal{O}(\frac{k^2}{\log N})O(kN)=O(logNk2)，对于k∼10000k\sim 10000k∼10000部分而言可能会比较危险，有TLE的可能。
怎么样优化掉这个O(nk)\mathcal{O}(nk)O(nk)的累加和？其实前面的Qingyu大佬有讲过，我们可以把所有的这些分式通分，变成：
eN=∑k=0N1k!=1N!∑k=0NANk=1N!∑k=0N∏t=N−k+1Nt\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 eN=k=0∑Nk!1=N!1k=0∑NANk=N!1k=0∑Nt=N−k+1∏Nt
这样，我们就可以分别计算分子和分母，再通过高精度除法得到最后的结果，这一过程中仅有最后的高精度除法是稍有误差的。前面使用高精度整数加法和与单精度的乘法，其时间复杂度大致与前面的以一致，为较小常数的O(k2logN)\mathcal{O}(\frac{k^2}{\log N})O(logNk2)；后面如果采取高精度除法而非模拟试除法，则可以优化到O(klogk)\mathcal{O}(k\log k)O(klogk)。经实验该方法本题可过。
本方法部分代码：
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罕见的高精度又带e\mathrm{e}e，那自然就不会放过这个机会测试函数exp\expexp啦。我们可以考虑实现exp(x)=ex\exp(x)=\mathrm{e}^xexp(x)=ex，再进行计算exp(1)\exp(1)exp(1)即可。那么我们接下来将着手分析并实现exp\expexp。
根据泰勒展开，我们可以将函数ex\mathrm{e}^xex展开成幂级数形式，并使用拉格朗日余项，我们有：
ex=∑k=0∞xkk!=∑k=0Nxkk!+RN(x)∣RN(x)∣=eθxxN+1(N+1)!≤ex(N+1)!\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)!} ex=k=0∑∞k!xk=k=0∑Nk!xk+RN(x)∣RN(x)∣=(N+1)!eθxxN+1≤(N+1)!ex
其中θ∈(0,1)\theta\in(0,1)θ∈(0,1)。由此得知如此展开以后的相对误差也不会超过1N!\frac{1}{N!}N!1。
同上面一样，虽然看起来本题可以计算了，然而仔细分析一下，高精度浮点数的乘法次数至少是O(N)\mathcal{O}(N)O(N)的，而高精度乘法的时间复杂度是O(klogk)\mathcal{O}(k\log k)O(klogk)的，放一起估计一下也有O(k2logk)\mathcal{O}(k^2\log k)O(k2logk)，那肯定会超时。虽然针对本题，后面会讲到如何优化掉这部分乘法，但总的来说还是要尽可能的降低这部分乘法的次数。
虽然与本题无关，但是这里还是简单的提一句：对于整个实数范围内的exp\expexp求值，我们可以将它规约到[0,1][0,1][0,1]的范围内。对于x<0x<0x<0的，有ex=1/e−x\mathrm{e}^{x}=1/\mathrm{e}^{-x}ex=1/e−x可以化为x>0x>0x>0的求值；对于x>1x>1x>1的，由常见的浮点数存储格式，可以认为x=x0×10Ex=x_0\times 10^Ex=x0×10E（其中x0∈(0,1)x_0\in (0,1)x0∈(0,1)并使用高精度整数来存储）(P.S. 这很像但并不是科学记数法，科学计数法可以规约到这种形式)，那么有exp(x)=exp(x0)10E\exp(x)=\exp(x_0)^{10^E}exp(x)=exp(x0)10E，可以化为(0,1)(0,1)(0,1)范围内的求值再做个高精度浮点数的正整数次幂得到。
参照PPT，我们可以进一步的将xxx的范围从(0,1](0,1](0,1]放缩到(0,10⌊−k⌋)(0,10^{\lfloor-\sqrt{k}\rfloor})(0,10⌊−k⌋)，然后仿照上面讨论x>1x>1x>1的部分那样提出10x10^{\sqrt{x}}10x放到exp\expexp的外面再做快速幂。对于x∈(0,10⌊−k⌋)x\in(0,10^{\lfloor-\sqrt{k}\rfloor})x∈(0,10⌊−k⌋)，有：
∣Rn(x)∣=eθxxn+1(n+1)!≤ex10−(n+1)⌊k⌋δ=Rn(x)ex<10−(n+1)⌊k⌋\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}}∣Rn(x)∣=(n+1)!eθxxn+1≤ex10−(n+1)⌊k⌋δ=exRn(x)<10−(n+1)⌊k⌋
则只需要取n=O(k)n=\mathcal{O}(\sqrt{k})n=O(k)，无论是exp\expexp的里面还是外面都只需要O(k)\mathcal{O}(\sqrt{k})O(k)次长度为kkk的高精度乘法，因此本方法的总的时间复杂度为O(k1.5logk)\mathcal{O}(k^{1.5}\log k)O(k1.5logk)。
虽然时间复杂度理论上满足了，然而实践操作中，如果是万进压位的高精度浮点数，对于k=10000k=10000k=10000的点，需要做大概500次的长度为4096~8192位的高精度乘法，这部分对常数的要求很大。我们可以采取一些比较通用的办法来进行优化。比如比较的经典的包括预处理循环卷积的单位根；又或者是对为循环卷积弄个缓存并适当调整乘法的位置，以保证不总是重复对同一个数组做DFT。当然也可以调整n=αkn=\alpha\sqrt{k}n=αk中的常数α\alphaα，平衡放缩前与放缩后的乘法次数。经过这些优化以后，我针对本题可以做到平均1.2s一个点，也许可以有更快的常数。
然而，针对本题，有一种更好的优化方法。如果你采用的是形如x=x0×10Ex=x_0\times 10^Ex=x0×10E这样的标准的存储法，你会发现无论是放缩前还是放缩后，xxx中的有效数位其实都只有一位。原因是不言自明的。因此在累加累乘的时候，只要在进行乘法时先判断两边的有效数位，其中的一边大概小于另一边的对数的级别的时候，就可以直接用朴素的按位数乘法累乘而不用经过循环卷积。这样可以减少至少一半的循环卷积次数，大大的提高运行效率。
部分核心代码如下：
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组成，其表示的数为base×10000pow\mathrm{base}\times10000^{\mathrm{pow}}base×10000pow。
扫一扫
245
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
