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) 式，也就是题目给出的公式计算 \mathrm{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 则需要 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.