P1729 计算e 的分治解法
TobyFlenderson
·
2024-06-06 22:54:30
·
题解
计算自然常数 \mathrm{e} 有多种办法，前人有很多研究。其中使用Maclaurin公式把 f(x)=\mathrm{e}^x 展开，并令 x=1 可以得到 \mathrm{e} 的无穷逼近式：
\mathrm{e}=1+\frac{1}{1!}+\frac{1}{2!}+\frac{1}{3!}+\cdots +\frac{1}{n!}+\cdots =1+\sum_{n=1}^{\infty}{\frac{1}{n!}}
记上式为 (1) 式。已有的研究表明 (1) 式收敛得很快[1] ，于是此题我们用 (1) 式，也就是题目给出的公式计算 \mathrm{e}。
拿到 (1) 式，首先要考虑的是 n 到底要取到多大才能满足题目所需的精度？这个问题 TBB_Nozomi 在他的题解中给出了详细的分析，这里引用他的结论：
要想使用 (1) 式计算 \mathrm{e} 的 n 位有效数字， n 应满足
10^n 两边取对数，有 n\ln 10<\sum_{k=1}^n{\ln k} 于是 n 可以编程求解，代码如下： int get_n(int d) { double f = d * log(10); double r = 0; int i = 0; do { i++; r += log(i); } while (r < f); return i; } 由于高精度除法的消耗很多，可以预测的是，如果直接用 (1) 式暴力计算 \mathrm{e} ，是一定会 TLE 的。于是我们要把 (1) 式通分： \mathrm{e}_n=1+\sum_{k=1}^n{\frac{1}{k!}}=1+\frac{1}{n!}\sum_{k=0}^{n-1}{A_{n}^{k}}， 使得原本要调用 n 次的高精度除法降为仅需 1 次。代码如下： void euler_tongfen(int d) { int n = get_n(d); mpz_class p = 1, q = 1; //mpz_class表示高精度整数类 for (int i = n; i > 1; i--) { q *= i; //高精度×单精度，时间复杂度为O(n) p += q; } mpf_class r = p; //mpq_class表示高精度浮点数类 r /= q; //高精度除法 ++r; //别忘了+1 } 显然上述算法的时间复杂度为 O(n^2) 。按照 TBB_Nozomi 的说法，只要高精度算法实现得没问题，上面的代码再加上一些格式化输出就能通过这题。下面是我用 TBB_Nozomi 的高精度整数板子，在他的题解代码的基础上修改得来的 AC 代码。 int main() { int k; cin >> k; tbb::LInt S = 1, P = 1; int N = get_n(k); for (int i = N; i > 0; i--) { P *= i; S += P; } S <<= k / 4 + 2; S /= P; string ans = S.print_str(); const char* out = ans.c_str(); putchar('2'); putchar('.'); putchar('\n'); for (int T = 1; T <= k; ++T) { putchar(out[T]); if (T % 50 == 0) putchar('\n'); else if (T % 10 == 0) putchar(' '); } return 0; } 这个链接是我的评测详情，开启了 O2 优化，总用时 255 ms。不开启 O2 则需要 556ms。 上述代码还有优化的可能： ①对于“高精度×单精度”的算法，它是 O(n) 的，在输入较小时它是十分快的，于是我们的优化思路是采用二分法对表达式分块处理，在表达式长度缩短到某一阈值的时候就采用通分的方式进行计算，然后使用时间复杂度低于 O(n^2) 的大整数乘法（FFT，Karatsuba，Toom Cook 3-Way）向上合并，最终得到结果。 ②对高精度除法进行优化，这就是另外一个问题了，可以查看P5432 A/B problem的题解。 如何进行分块处理呢？为了方便讨论，记 E\left( n,m \right) =\frac{1}{n}+\frac{1}{n\left( n+1 \right)}+\cdots +\frac{1}{n\left( n+1 \right) \left( n+2 \right) \cdots m} 于是 \mathrm{e}=1+E(1,+\infty) 。令 r=\lfloor \frac{n+m}{2} \rfloor，那么对于 E(n,m) 有 E\left( n,m \right) =E\left( n,r \right) +\frac{1}{n\left( n+1 \right) \cdots r}\cdot E\left( r+1,m \right) 再令 \frac{p_1}{q_1}=E\left( n,r \right)，\frac{p_2}{q_2}=E\left( r+1,m \right)，显然 q_1=n(n+1)\cdots r，于是 E\left( n,m \right) =\frac{p_1}{q_1}+\frac{1}{q_1}\cdot \frac{p_2}{q_2}=\frac{p_1q_2+p_2}{q_1q_2} 算法如下图所示 每做一次合并，需要两次乘法，一次加法。显然，若“高精×高精”的算法复杂度仍为 O(n^2) ，那么总体的复杂为 T\left( n \right) =\frac{n^2}{2}+\frac{n^2}{4}+\frac{n^2}{8}+\cdots =n^2\left( \frac{1}{2}+\frac{1}{4}+\frac{1}{8}+\cdots \right) =O\left( n^2 \right) 与原算法基本持平，甚至还要稍慢一些。于是只要“高精×高精“的算法复杂度低于 O(n^2) ，那么总体的复杂度就比原来的更优。 那么只剩下最后一个问题了，阈值的选取，也就是 E(n,m) 多短我们才开始用通分法计算呢？遗憾的是，不同的处理器，不同的算法实现效率都会影响阈值的选取，要找到阈值最简单的方法是通过暴力搜索的方式得到。在我的机器上这个值大约在 120 左右。在洛谷上使用阈值为 128 的新算法用时 161ms，评测详情。 以下是新的 AC 代码。 int MIN_SPLIT = 128; static void euler_split(int n, int m, LInt& p, LInt& q) { if (m - n < MIN_SPLIT) { p = 1; q = 1; for (int i = m; i > n; i--) { q *= i; p += q; } q *= n; return; } LInt p1, p2, q1, q2; euler_split(n, (n + m) >> 1, p1, q1); euler_split((n + m + 2) >> 1, m, p2, q2); p = p1 * q2 + p2; q = q1 * q2; } int main() { int k; cin >> k; LInt p, q; int n = get_n(k); euler_split(1, n, p, q); p += q; p <<= k / 4 + 2; p /= q; string ans = p.print_str(); const char* out = ans.c_str(); putchar('2'); putchar('.'); putchar('\n'); for (int T = 1; T <= k; ++T) { putchar(out[T]); if (T % 50 == 0) putchar('\n'); else if (T % 10 == 0) putchar(' '); } return 0; } [1] 朱德凯, 苏岐芳. 多种计算方法下自然常数e的误差和收敛速度比较[J]. 台州学院学报, 2022, 44(03): 11-16. DOI: 10.13853/j.cnki.issn. 1672-3708. 2022. 03. 003.