P4726 【模版】多项式指数函数（多项式exp）

P4726 【模版】多项式指数函数（多项式exp）

题面

题目链接

解法

定义：给定多项式$A(x)$，求$B(x)$使得$B(x) \equiv e^{A(x)} \pmod{x^n}$。

根据题目有$B(x) = e^{A(x)}$，两边同时求导：

$$\begin{aligned} B'(x) =& e^{A(x)} \times A'(x) \\ B'(x) =& B(x) \times A'(x) \\ B(x) =& \int B(x) \times A'(x) \end{aligned}$$

使用分治NTT暴力搞，$O(n \log^2 n)$（代码不贴了，毕竟不是最优）。

牛顿迭代法

设想我们知道一个数$a$，要求$\sqrt{a}$精确到整数。

首先能想到二分，但是当$a$非常大，需要二分的次数就需要$\log_2 n$次。

于是使用牛顿迭代法求解。设$f(x)=x^2-a$，要求的就是$f(x)$的零点。

假设我们求出了一个近似值$x_0$，就可以过$(x_0,f(x_0))$作$f(x)$的切线，取切线与$x$轴的交点作为新的$x_0$。这样需要的次数是远小于直接二分的。

具体地，当我们有一个近似值$x_0$，其关于$f(x)$的切线方程如下：

$$y=f'(x_0)(x-x_0)+f(x_0)$$

求出切线与$x$轴的交点$x$：

$$x=x_0-\frac{f(x_0)}{f'(x_0)}$$

以上是二维坐标系中的牛顿迭代，而多项式可以看成是多维坐标系中的点，总之牛顿迭代仍然使用。

想要求一个函数$G(x)$使得$F(G(x))\equiv 0$。

使用上面的式子，我们呢每次令

$$G(x)=G_0(x)-\frac{F(G_0(x))}{F'(G_0(x))}$$

经过严谨的证明（我不会证），每次迭代可以使精度翻倍，即当$F(G_0(x))\equiv 0 \pmod{x^\frac{n}{2}}$，迭代后$F(G(x))\equiv 0 \pmod{x^n}$。

带到这道题中，我们首先由

$$\begin{aligned} B(x) \equiv e^{A(x)} &\pmod{x^n} \\ \ln B(x) \equiv A(x) &\pmod{x^n} \\ \ln B(x) - A(x) \equiv 0 &\pmod{x^n} \\ \end{aligned}$$

我们让$F(G(x))=\ln G(x) - A(x) \equiv 0$即可。

将$F(G(x))$求导，得到：

$$F'(G(x)) = \frac{1}{G(x)}$$

因为对于函数$F(G(x))$而言，自变量是$G(x)$而不是$x$，其中的$A(x)$相当于是常数，因而直接被忽略。

将导数带进上式得到：

$$\begin{aligned} G(x) =& G_0(x) - \frac{\ln G_0(x) - A(x)}{\frac{1}{G_0(x)}} \\ =& G_0(x)(1 + A(x) - \ln G_0(x)) \end{aligned}$$

递归计算，边界条件是$[x^0]B(x)=1$。

所以时间复杂度为

$$T(n) = T(\frac{n}{2}) + O(n \log n) = O(n \log n)$$

AC代码

/**
 * @file:           P4726.cpp
 * @author:         yaoxi-std
 * @url:            https://www.luogu.com.cn/problem/P4726
*/
// #pragma GCC optimize ("O2")
// #pragma GCC optimize ("Ofast", "inline", "-ffast-math")
// #pragma GCC target ("avx,sse2,sse3,sse4,mmx")
#include <bits/stdc++.h>
using namespace std;
#define int long long
#define resetIO(x) \
    freopen(#x ".in", "r", stdin), freopen(#x ".out", "w", stdout)
#define debug(fmt, ...) \
    fprintf(stderr, "[%s:%d] " fmt "\n", __FILE__, __LINE__, ##__VA_ARGS__)
template <class _Tp>
inline _Tp& read(_Tp& x) {
    bool sign = false;
    char ch = getchar();
    long double tmp = 1;
    for (; !isdigit(ch); ch = getchar())
        sign |= (ch == '-');
    for (x = 0; isdigit(ch); ch = getchar())
        x = x * 10 + (ch ^ 48);
    if (ch == '.')
        for (ch = getchar(); isdigit(ch); ch = getchar())
            tmp /= 10.0, x += tmp * (ch ^ 48);
    return sign ? (x = -x) : x;
}
template <class _Tp>
inline void write(_Tp x) {
    if (x < 0)
        putchar('-'), x = -x;
    if (x > 9)
        write(x / 10);
    putchar((x % 10) ^ 48);
}
const int MAXN = 1 << 19;
const int INF = 0x3f3f3f3f3f3f3f3f;
const int MOD = 998244353;
namespace maths {
int add(int x, int y) {
    x += y;
    return x >= MOD ? x - MOD : x;
}
int sub(int x, int y) {
    x -= y;
    return x < 0 ? x + MOD : x;
}
int qpow(int x, int y, int p = MOD) {
    int ret = 1;
    for (; y; y >>= 1, x = x * x % p)
        if (y & 1)
            ret = ret * x % p;
    return ret;
}
template <class _Tp>
void change(_Tp* f, int len) {
    static int rev[MAXN] = {};
    for (int i = 0; i < len; ++i) {
        rev[i] = rev[i >> 1] >> 1;
        if (i & 1)
            rev[i] |= len >> 1;
    }
    for (int i = 0; i < len; ++i)
        if (i < rev[i])
            swap(f[i], f[rev[i]]);
}
void ntt(int* f, int len, int on) {
    change(f, len);
    for (int h = 2; h <= len; h <<= 1) {
        int gn = qpow(3, (MOD - 1) / h);
        for (int j = 0; j < len; j += h) {
            int g = 1;
            for (int k = j; k < j + h / 2; ++k) {
                int u = f[k], t = g * f[k + h / 2] % MOD;
                f[k] = add(u, t), f[k + h / 2] = sub(u, t);
                g = g * gn % MOD;
            }
        }
    }
    if (on == -1) {
        reverse(f + 1, f + len);
        int inv = qpow(len, MOD - 2);
        for (int i = 0; i < len; ++i)
            f[i] = f[i] * inv % MOD;
    }
}
void poly_multiply(const int* f, int n, const int* g, int m, int* ans) {
    static int tf[MAXN] = {}, tg[MAXN] = {};
    int len = 1;
    while (len < n + m - 1)
        len <<= 1;
    fill(ans, ans + len, 0);
    copy(f, f + n, tf);
    fill(tf + n, tf + len, 0);
    copy(g, g + m, tg);
    fill(tg + m, tg + len, 0);
    ntt(tf, len, 1);
    ntt(tg, len, 1);
    for (int i = 0; i < len; ++i)
        tf[i] = tf[i] * tg[i] % MOD;
    ntt(tf, len, -1);
    copy(tf, tf + n + m - 1, ans);
}
void poly_inverse(const int* f, int n, int* ans) {
    static int tmp[MAXN] = {};
    int len = 1;
    while (len < n)
        len <<= 1;
    fill(ans, ans + len + len, 0);
    ans[0] = qpow(f[0], MOD - 2);
    for (int h = 2; h <= len; h <<= 1) {
        copy(f, f + h, tmp);
        fill(tmp + h, tmp + h + h, 0);
        ntt(tmp, h + h, 1);
        ntt(ans, h + h, 1);
        for (int i = 0; i < h + h; ++i)
            ans[i] = ans[i] * (2 - ans[i] * tmp[i] % MOD + MOD) % MOD;
        ntt(ans, h + h, -1);
        fill(ans + h, ans + h + h, 0);
    }
}
void poly_derivation(const int* f, int n, int* ans) {
    for (int i = 0; i < n - 1; ++i)
        ans[i] = f[i + 1] * (i + 1) % MOD;
    ans[n - 1] = 0;
}
void poly_integral(const int* f, int n, int* ans) {
    for (int i = n; i >= 1; --i)
        ans[i] = f[i - 1] * qpow(i, MOD - 2) % MOD;
    ans[0] = 0;
}
void poly_ln(const int* f, int n, int* ans) {
    static int tf[MAXN] = {}, tg[MAXN] = {};
    poly_derivation(f, n, tf);
    poly_inverse(f, n, tg);
    poly_multiply(tf, n - 1, tg, n, ans);
    poly_integral(ans, n - 1, ans);
    fill(ans + n, ans + n + n, 0);
}
void poly_exp(const int* f, int n, int* ans) {
    static int tf[MAXN] = {}, tg[MAXN] = {};
    ans[0] = 1, ans[1] = 0;
    for (int h = 2; h <= (n << 1); h <<= 1) {
        poly_ln(ans, h, tf);
        for (int i = 0; i < h; ++i)
            tf[i] = add((i == 0), sub(f[i], tf[i]));
        copy(ans, ans + h, tg);
        poly_multiply(tf, h, tg, h, ans);
    }
}
};  // namespace maths
using namespace maths;
int n, a[MAXN], b[MAXN], c[MAXN];
signed main() {
    read(n);
    for (int i = 0; i < n; ++i)
        read(a[i]);
    poly_exp(a, n, b);
    for (int i = 0; i < n; ++i)
        write(b[i]), putchar(" \n"[i == n - 1]);
    return 0;
}

此后，我决定一改封装多项式的坏习惯，转而使用面向过程的函数式编程，以防止多项式传递时额外的内存开销。

同时，近一个月的多项式学习，让我深刻认识到自己数学知识的不足。因此，这个寒假需要做的事不言而喻。