Luogu 4726 【模板】多项式指数函数

补补补……

这个题的解法让我认识到了泰勒展开的美妙之处。

泰勒展开

泰勒展开就是用一个多项式型的函数去逼近一个难以准确描述的函数。

有公式

$$f(x)approx g(x) = g(x_0) + frac{g'(x_0)}{1!}(x - x_0) + frac{g^{(2)}(x_0)}{2!}(x - x_0)^2 + cdots + frac{g^{(n)}(x_0)}{n!}(x - x_0)^n$$

在这里$g^{(n)}$表示$g$的$n$阶导。

在$0$这个点的泰勒展开$(x_0 = 0)$叫做麦克劳林级数，利用这个东西可以很精确地去逼近原来的函数。

比如，

$$f(x) = e^x approx 1 + frac{x}{1!} + frac{x^2}{2!} + frac{x^3}{3!} + frac{x^4}{4!} + cdots = sum_{i = 0}^{infty}frac{x^i}{i!}$$

$$f(x) = ln(1 - x) approx 0 - frac{x}{1} - frac{x^2}{2} - frac{x^3}{3} - cdots = -sum_{i = 1}^{infty}frac{x^i}{i}$$

一些定义

有了泰勒展开之后我们可以定义一些看上去很难定义的东西，比如这个多项式取$ln$:

$$ln(1 - A(x)) = -sum_{i = 1}^{infty}frac{A(x)^i}{i}$$

这也解释了为什么给出的多项式常数项一定是$1$。

同理，多项式$exp$的定义是这样子的:

$$exp(A(x)) = sum_{i = 0}^{infty}frac{A(x)^i}{i!}$$

可以发现这样子定义之后就好理解很多了。

牛顿迭代

牛顿迭代是解决多项式问题的利器，可以一口气解决好多推式子问题。

牛顿法求解方程的近似解在高中课本里面有提到过，它也可以用来解决多项式的问题。

现在我们尝试用另一种手段得到它。

要求$G(F(x)) equiv 0 (mod x^n)$，通过题目给出的特定的$G(x)$，我们可以算出$F(x)$。

首先我们尝试解决常数项的问题，一般来说，这个问题都非常简单，比如本题中可以直接求出常数项为$1$。

现在假设已经求出了$F_0(x)$满足$G(F_0(x)) equiv 0(mod x^{left lceil frac{n}{2} ight ceil})$，我们尝试求$F(x)$满足$G(F(x)) equiv 0 (mod x^n)$。

可以对$G(x)$在$F_0(x)$处进行泰勒展开，得到

$$G(F(x)) = G(F_0(x)) + G'(F_0(x))(F(x) - F_0(x)) + G''(F_0(x))(F(x) - F_0(x))^2 + cdots$$

这是在$(mod x^n)$次意义下进行的，在多项式求逆的时候已经证明过了$(F(x) - F_0(x))^n$在$n geq 2$的时候模$x^n$为$0$，后面就可以不用写下去了。

注意到$G(F(x)) equiv 0 (mod x^n)$，有

$$F(x) equiv F_0(x) - frac{G(F_0(x))}{G'(F_0(x))}(mod x^n)$$

这个式子就是牛顿迭代的式子了。

鼓掌~~~

考虑一下在这个题中怎么取这个$G(x)$。

题目要求

$$B(x) equiv e^{A(x)} (mod x^n)$$

发现这个$e^x$并不是很方便，所以两边取一下$ln$，移项，

$$ln(B(x)) - A(x) equiv 0 (mod x^n)$$

那就记$G(F(x)) = lnF(x) - A(x)$。

因为$A(x)$和$F(x)$在这里等价于两个数，$F(x)$是自变量，所以$A(x)$就看成常数，在求导的时候可以直接消掉。

整理一下就得到了多项式$exp$的式子:

$$F(x) = F_0(x)(1 - lnF_0(x) + A(x))$$

左转去复制各种模板过来，于是就可以愉快地递归求解了！

时间复杂度仍然是$O(nlogn)$。

在吉老师博客里看到了关于这堆东西常数的吐槽，感觉妙不可言。

Code：

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;

const int N = 1 << 18;

int n;
ll f[N], g[N];

namespace Poly {
    const int L = 1 << 18;
    const ll P = 998244353LL;

    int lim, pos[L];
    ll f[L], g[L], h[L], tmp[L];
    
    inline ll fpow(ll x, ll y) {
        ll res = 1;
        for (x %= P; y > 0; y >>= 1) {
            if (y & 1) res = res * x % P;
            x = x * x % P;
        }
        return res;
    }
    
    inline void prework(int len) {
        int l = 0;
        for (lim = 1; lim < len; lim <<= 1, ++l);
        for (int i = 0; i < lim; i++)
            pos[i] = (pos[i >> 1] >> 1) | ((i & 1) << (l - 1));
    }
    
    inline void ntt(ll *c, int opt) {
        for (int i = 0; i < lim; i++)
            if (i < pos[i]) swap(c[i], c[pos[i]]);
        for (int i = 1; i < lim; i <<= 1) {
            ll wn = fpow(3, (P - 1) / (i << 1));
            if (opt == -1) wn = fpow(wn, P - 2);
            for (int len = i << 1, j = 0; j < lim; j += len) {
                ll w = 1;
                for (int k = 0; k < i; k++, w = w * wn % P) {
                    ll x = c[j + k], y = w * c[j + k + i] % P;
                    c[j + k] = (x + y) % P, c[j + k + i] = (x - y + P) % P;
                }
            }
        }
        
        if (opt == -1) {
            ll inv = fpow(lim, P - 2);
            for (int i = 0; i < lim; i++) c[i] = c[i] * inv % P;
        }
    }
    
    void inv(ll *a, ll *b, int len) {
        if (len == 1) {
            b[0] = fpow(a[0], P - 2);
            return;
        }
        
        inv(a, b, (len + 1) >> 1);
        prework(len << 1);
        for (int i = 0; i < lim; i++) f[i] = g[i] = 0;
        for (int i = 0; i < len; i++) f[i] = a[i], g[i] = b[i];
        ntt(f, 1), ntt(g, 1);
        for (int i = 0; i < lim; i++)
            g[i] = g[i] * (2LL - g[i] * f[i] % P + P) % P;
        ntt(g, -1);
        
        for (int i = 0; i < len; i++) b[i] = g[i];
    }
    
    inline void direv(ll *c, int len) {
        for (int i = 0; i < len - 1; i++) c[i] = c[i + 1] * (i + 1) % P;
        c[len - 1] = 0;
    }
    
    inline void integ(ll *c, int len) {
        for (int i = len - 1; i > 0; i--) c[i] = c[i - 1] * fpow(i, P - 2) % P;
        c[0] = 0;
    }
    
    inline void ln(ll *a, ll *b, int len) {
        for (int i = 0; i < len; i++) h[i] = 0;
        inv(a, h, len);
                
        for (int i = 0; i < len; i++) b[i] = a[i];
        direv(b, len);    

        prework(len << 1);
        ntt(h, 1), ntt(b, 1);
        for (int i = 0; i < lim; i++) b[i] = b[i] * h[i] % P;
        ntt(b, -1);
        
        integ(b, len);
//        for (int i = 0; i < lim; i++) h[i] = 0;
    }
    
    ll F[L], G[L];
    void exp(ll *a, ll *b, int len) {
        if (len == 1) {
            b[0] = 1;
            return;
        }
        exp(a, b, (len + 1) >> 1);
        
        ln(b, F, len);
        F[0] = (a[0] % P + 1 - F[0] + P) % P;
        for (int i = 1; i< len; i++) F[i] = (a[i] - F[i] + P) % P;
        
        prework(len << 1);
        for (int i = len; i < lim; i++) F[i] = 0;
        for (int i = 0; i < lim; i++) G[i] = 0;
        for (int i = 0; i < len; i++) G[i] = b[i];
        ntt(F, 1), ntt(G, 1);
        for (int i = 0; i < lim; i++) F[i] = F[i] * G[i] % P;
        ntt(F, -1);
        
        for (int i = 0; i < lim; i++) b[i] = F[i];
    }
    
};

template <typename T>
inline void read(T &X) {
    X = 0; char ch = 0; T op = 1;
    for (; ch > '9'|| ch < '0'; ch = getchar())
        if (ch == '-') op = -1;
    for (; ch >= '0' && ch <= '9'; ch = getchar())
        X = (X << 3) + (X << 1) + ch - 48;
    X *= op;
}

int main() {
    read(n);
    for (int i = 0; i < n; i++) read(f[i]);
    Poly :: exp(f, g, n);
    for (int i = 0; i < n; i++)
        printf("%lld%c", g[i], " 
"[i == n - 1]);
    return 0;
}

 一点点推广

取不同的$G(x)$可以得到不同的结果，简单推导出答案的式子。

比如：

多项式求逆，

$$G(F(x)) = F(x)A(x) - 1$$

多项式开方，

$$G(F(x)) = F(x)^2 - A(x)$$

……

还有一些应该会在做题的时候补全。