洛谷 P1729 计算e 的分治解法

计算自然常数 $\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)$ 式，也就是题目给出的公式计算e。

拿到 $(1)$ 式，首先要考虑的是 $n$ 到底要取到多大才能满足题目所需的精度？这个问题 TBB_Nozomi 在他的题解中给出了详细的分析，这里引用他的结论：

要想使用 $(1)$ 式计算 $\mathrm{e}$ 的 $n$ 位有效数字， $n$ 应满足
$$10^n<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则需要 $556$ms。

上述代码还有优化的可能：

①对于“高精度×单精度”的算法，它是 $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_1{q}_2+p_2}{q_1{q}_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$ 的新算法用时 $161$ms，评测详情。

以下是新的 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.