[干饭现场] 多项式全家桶

-1. 无关紧要的话

\(\text{kfc}\) 全家桶都只够两个人吃，我这个全家桶不收录所有模板也不过分吧😶‍🌫️。

顺便推销一波自己的 \(\mathtt{FFT}\) 和 \(\mathtt{NTT}\).

0. 泰勒展开与多项式牛顿迭代

\(\color{red}{\text{Disclaimer}}\)：由于辣鸡博主数学很烂，所以可能会有很多不严谨的地方，如果哪些地方有问题请一定教教我 \(\text{sto}\).

0.1. 泰勒展开

0.1.1. 展开方式

泰勒展开是将一种在 \(x=x_0\) 处具有 \(n\) 阶导数的函数 \(f(x)\) 利用关于 \((x-x_0)\) 的 \(n\) 次来逼近函数多项式的方法。

若函数 \(f(x)\) 在包含 \(x_0\) 的某个闭区间 \([a,b]\) 上具有 \(n\) 阶导数，且在开区间 \((a,b)\) 上具有 \((n+1)\) 阶导数，则对 \([a,b]\) 上任意一点 \(x\)，成立下式

\[\begin{aligned} f(x)=\sum _{i=0}^{n}\frac{f^{(i)}(x_0)}{i!}\cdot (x-x_0)^i\end{aligned}+R_n(x) \]

等号后的多项式称为函数 \(f(x)\) 在 \(x_0\) 处的泰勒展开式，剩余的 \(R_n(x)\) 是泰勒公式的余项，是 \((x-x_0)^n\) 的高阶无穷小，相当于泰勒展开式对原函数 \(f(x)\) 的一些相对而言极小的误差。每一次展开都是一次逼近。

当然，如果函数可以无限次求导的话，泰勒展开也可以无限次展开。但事实上对于我们所研究的有限次多项式，经过有限次求导就会变成 \(0\)，展开次数仍然是有限的。

0.1.2. 一些常见的展开

\(\sin x\) 的 \(0,1,2,3\) 阶导数为 \(\sin x,\cos x,-\sin x,-\cos x\)。这是一个循环。

取 \(x_0=0\)，泰勒展开有

\[\begin{align} \sin x&=\sum_{i=0}^{\infty}\frac{\sin^{(i)}(0)}{i!}\cdot x^i\\ &=\sum_{k=0}^{\infty}\frac{\sin^{(2k)}(0)}{2k!}\cdot x^{2k}+\sum_{k=0}^{\infty}\frac{\sin^{(2k+1)}(0)}{(2k+1)!}\cdot x^{2k+1}\\ &=\sum_{k=0}^{\infty}\frac{(-1)^k}{(2k+1)!}\cdot x^{2k+1} \end{align} \]
\(\ln (1+x)= \displaystyle \sum_{n\geqslant 1}(-1)^{n+1}\cdot \dfrac{x^n}{n}\)，同时可以发现 \((1+x)\) 可以很好地消去 \(n=0\) 的情况。

同样地，\(\ln (1-x)=\ln (1+(-x))=- \displaystyle \sum_{n\geqslant 1}\dfrac{x^n}{n}\).
其实是 \(\ln(1+x)\) 展开的拓展，把系数移了一位，\(\dfrac{1}{1+x}= \displaystyle \sum_{n\geqslant 0} (-1)^nx^n\).
还可以用泰勒展开证明二项式定理：\((1+x)^a= \displaystyle \sum_{n\geqslant 0}\text{C}(a,n)\cdot x^n\).
另外还有：\(\dfrac{1}{(1-x)^{a+1}}= \displaystyle \sum_{n\geqslant 0}\dbinom{a+n}{n}\cdot x^n\)，特别地，当 \(a=1\) 时存在 \(\dfrac{1}{(1-x)^2}= \displaystyle \sum_{n\geqslant 0}(n+1)\cdot x^n\).
\(\text{e}^x= \displaystyle \sum_{n\geqslant 0} \dfrac{x^n}{n!}\)，另外 \({\rm e}^{px}= \displaystyle \sum_{n\geqslant 0}\dfrac{p^nx^n}{n!}\).

0.2. 多项式牛顿迭代

0.2.1. 问题描述

给定多项式 \(g(x)\)，已知有 \(f(x)\) 满足

\[g(f(x))\equiv 0\pmod{x^n} \]

求出模 \(x^n\) 意义下的 \(f(x)\).

其中 \(g(f(x))=\sum_{i=0}^{\infty}g_i\cdot f^i(x)\)。需要特别注意的是，\(g\) 是关于 \(f(x)\) 的函数，而不是关于 \(x\) 的函数。

0.2.2. 牛顿迭代

考虑倍增。首先当 \(n=1\) 时，需要预先算出对应 \(f(x)\)。记 \(f_0(x)\) 为在模 \(x^{n/2}\) 意义下的解，那么有

\[\begin{align} g(f_0(x))&\equiv 0\pmod{x^{n/2}}\\ g(f(x))&\equiv 0\pmod{x^{n/2}} \end{align} \]

移项得

\[f(x)-f_0(x)\equiv 0\pmod{x^{n/2}} \]

后面有类似的证明，这里不再赘述，反正可以得到

\[(f(x)-f_0(x))^2\equiv 0\pmod{x^n} \]

\(\rm ok\) 接下来将 \(g(f(x))\) 在 \(f_0(x)\) 处泰勒展开^[1]

\[g(f(x))=\sum_{i=0}^{\infty}\frac{g^{(i)}(f_0(x))}{i!}\cdot (f(x)-f_0(x))^i \]

由上文可知，\(i\ge 2\) 的项都可以被舍去^[2]，解一下方程即可得到

\[f(x)\equiv f_0(x)-\frac{g(f_0(x))}{g'(f_0(x))}\pmod{x^n} \]

1. 多项式求逆

1.1. 问题描述

给定多项式 \(f\)，求多项式 \(g\) 满足 \(f*g\equiv 1\ (\text{mod}\ x^n)\).

1.2. 解法

先开始属实没看懂这是什么意思。其实翻译过来就是卷积后满足前 \(n\) 项系数除了常数项为 \(1\)，其余均为 \(0\).

可以考虑倍增求解。若已求得 \(g'\) 满足

\[f*g'\equiv 1\ (\text{mod}\ x^{\frac{n}{2}}) \]

显然要求的 \(g\) 也要满足

\[f*g\equiv 1\ (\text{mod}\ x^{\frac{n}{2}}) \]

故

\[\begin{align} f*g*g&\equiv f*g'*g\ (\text{mod}\ x^{\frac{n}{2}})\\ g&\equiv g'\ (\text{mod}\ x^{\frac{n}{2}})\\ g-g'&\equiv 0\ (\text{mod}\ x^{\frac{n}{2}}) \end{align} \]

将多项式 \(g-g'\) 与其本身做卷积，由于其 \(0\) 到 \(n/2-1\) 项系数都为 \(0\)，容易发现卷积完 \(0\) 到 \(n-1\) 项系数都为 \(0\)（考虑 \(n-1\) 最平均的分解也是 \((n/2-1)+(n/2)\) 了）。另外，\(n/2\) 需要 向上取整。

即

\[(g-g')^2\equiv 0\ (\text{mod}\ x^n) \]

拆项有

\[\begin{align} g^2-2gg'+(g')^2&\equiv 0\ (\text{mod}\ x^{n})\\ g-2g'+(g')^2*f&\equiv 0\ (\text{mod}\ x^{n})\\ g&\equiv 2g'-(g')^2*f \ (\text{mod}\ x^n) \end{align} \]

递归求解即可，复杂度 \(\mathcal O(n\log n)\)。这是因为里面层数加起来的复杂度也不及最外一层的复杂度。

上面就是朴素地推导多项式求逆的过程，接下来以其为例使用牛顿迭代再次推导。

假设 \(h(x)\) 是需要求逆的函数，\(f(x)\) 为它的逆，记 \(g(f(x))=\frac{1}{f(x)}-h(x)\)，容易发现这个式子在模 \(x^n\) 意义下是同余于 \(0\) 的。直接套用牛顿迭代得

\[\begin{align} f(x)&\equiv f_0(x)-\frac{\frac{1}{f_0(x)}-h(x)}{-\frac{1}{f_0^2(x)}}\pmod{x^n}\\ &\equiv 2f_0(x)-f_0^2(x)\cdot h(x)\pmod{x^n} \end{align} \]

这里再说一下为啥 \(g(f_0(x))\) 求导是 \(-\frac{1}{f_0^2(x)}\)，这是因为 \(g\) 是关于 \(f_0(x)\) 的函数！所以 \(h(x)\) 相当于常数项。

可以发现，我们只需要配凑一个关于原函数和目标函数同余零的式子，同时将目标函数设为 \(g\) 的自变量，就可以使用牛顿迭代求解多项式的各类运算。

1.3. 代码实现

先开始想着写递推版会快一些，结果跑出来发现比之前写的递归版慢了一倍😶‍🌫️。后来发现是因为递归版有个很大的优化 —— 从给定长度 \(n\) 开始除以二。而递推版是从 \(2\) 开始倍增，当长度不小于 \(n\) 时才开始卷，此时的上界就变成 \(8n\) 了！~~虽然改完还是要跑五百多毫秒~~。

另外要注意的代码细节是从 \(2\) 开始卷而不是从 \(1\) 开始卷，因为从 \(g'\) 递推到 \(g\) 是在模 \(x^n\) 意义下而不是在模 \(x^{n/2}\) 意义下完成的。

# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)

template <class T>
inline T read(const T sample) {
    T x=0; char s; bool f=0;
    while(!isdigit(s=getchar())) f|=(s=='-');
    for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
    return f? -x: x;
}
template <class T>
inline void write(T x) {
    static int writ[50], w_tp=0;
    if(x<0) putchar('-'), x=-x;
    do writ[++w_tp]=x-x/10*10, x/=10; while(x);
    while(putchar(writ[w_tp--]^48), w_tp);
}

# include <cstring>
# include <iostream>
using namespace std;

const int maxn = 4e5+5;
const int mod = 998244353, g = 3;

inline int inv(int x,int y=mod-2,int r=1) {
	for(; y; y>>=1, x=1ll*x*x%mod)
		if(y&1) r=1ll*r*x%mod;
	return r;
} const int ig = inv(g);
inline int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
inline int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }

int rev[maxn],n,F[maxn],G[maxn];

inline void ntt(int* f,int lim,bool opt=0) {
    static int w, wn, tmp;
    for(int i=0;i<lim;++i) 
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1; mid<lim; mid<<=1) {
        wn = inv(opt?ig:g,(mod-1)/(mid<<1)), w=1;
        for(int i=0; i<lim; i+=(mid<<1), w=1)
            for(int j=0; j<mid; ++j, w=1ll*w*wn%mod)
                tmp = 1ll*w*f[i|j|mid]%mod,
                f[i|j|mid] = dec(f[i|j],tmp),
                f[i|j] = inc(f[i|j],tmp);
    }
}

inline void polyInv(int* g,const int* h,int n) {	
    static int f[maxn], Inv, bit[30], b_len=0, lim, B, w; 
    while(n^1) bit[++ b_len] = n, n = n+1>>1; g[0] = inv(h[0]);
    for(; b_len; --b_len) {
        w=bit[b_len], B=0, lim=1;
        while(lim<(w<<1)) lim<<=1, ++B;
        memcpy(f,h,sizeof(int)*w);
        memset(f+w,0,sizeof(int)*(lim-w));
        for(int i=0;i<lim;++i) 
            rev[i] = (rev[i>>1]>>1)|((i&1)<<B-1);
        ntt(g,lim), ntt(f,lim);
        for(int i=0;i<lim;++i)
            g[i] = 1ll*g[i]*dec(2,1ll*g[i]*f[i]%mod)%mod;
        ntt(g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) g[i] = 1ll*g[i]*Inv%mod; 
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

int main() {
    n=read(9); 
    for(int i=0;i<n;++i) F[i]=read(9);
    polyInv(G,F,n);
    for(int i=0;i<n;++i) print(G[i],' '); puts("");
	return 0;
}

2. 多项式开根

2.1. 问题描述

给定多项式 \(f\)，求多项式 \(g\) 满足 \(g^2\equiv f\ (\text{mod}\ x^n)\).

2.2. 解法

可以考虑倍增求解。若已求得 \(g'\) 满足

\[(g')^2\equiv f\ (\text{mod}\ x^{\frac{n}{2}}) \]

显然要求的 \(g\) 也要满足

\[g^2\equiv f\ (\text{mod}\ x^{\frac{n}{2}}) \]

故（由平方差公式）

\[g-g'\equiv 0\ (\text{mod}\ x^{\frac{n}{2}}) \]

同上可得

\[(g-g')^2\equiv 0\ (\text{mod}\ x^n) \]

拆项有

\[g^2-2gg'+(g')^2\equiv 0\ (\text{mod}\ x^{n}) \]

\[f-2gg'+(g')^2\equiv 0\ (\text{mod}\ x^{n}) \]

\[g\equiv \frac{f+(g')^2}{2g'} \ (\text{mod}\ x^n) \]

复杂度同样也是 \(\mathcal O(n\log n)\)，因为求逆与开根的运算实际上是并行的。

2.3. 代码实现

# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)

template <class T>
inline T read(const T sample) {
    T x=0; char s; bool f=0;
    while(!isdigit(s=getchar())) f|=(s=='-');
    for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
    return f? -x: x;
}
template <class T>
inline void write(T x) {
    static int writ[50], w_tp=0;
    if(x<0) putchar('-'), x=-x;
    do writ[++w_tp]=x-x/10*10, x/=10; while(x);
    while(putchar(writ[w_tp--]^48), w_tp);
}

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

const int maxn = 4e5+5;
const int mod = 998244353, g = 3;

inline int inv(int x,int y=mod-2,int r=1) {
	for(; y; y>>=1, x=1ll*x*x%mod)
		if(y&1) r=1ll*r*x%mod;
	return r;
} const int ig = inv(g), iv2 = inv(2);
inline int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
inline int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }

int rev[maxn],n,F[maxn],G[maxn];

inline void ntt(int* f,int lim,bool opt=0) {
    static int tmp, w, wn;
    for(int i=0;i<lim;++i)
        if(i<rev[i]) swap(f[i],f[rev[i]]); 
    for(int mid=1; mid<lim; mid<<=1) {
        wn = inv(opt? ig: g, (mod-1)/(mid<<1)), w=1;
        for(int i=0; i<lim; i+=(mid<<1), w=1)
            for(int j=0; j<mid; ++j, w=1ll*w*wn%mod) 
                tmp = 1ll*w*f[i|j|mid]%mod,
                f[i|j|mid] = dec(f[i|j],tmp),
                f[i|j] = inc(f[i|j],tmp);
    }
}

inline void polyInv(int* g,int* h,int n) {
    static int f[maxn],bit[30],b_len=0,B,lim,w,Inv;
    while(n^1) bit[++b_len] = n, n = n+1>>1; g[0]=inv(h[0]);
    for(; b_len; --b_len) {
        w=bit[b_len], B=0, lim=1;
        while(lim<(w<<1)) lim<<=1, ++B;
        for(int i=0;i<lim;++i) 
            rev[i] = (rev[i>>1]>>1)|((i&1)<<B-1);
        memcpy(f,h,sizeof(int)*w);
        memset(f+w,0,sizeof(int)*(lim-w));
        ntt(f,lim), ntt(g,lim);
        for(int i=0;i<lim;++i)
            g[i] = 1ll*g[i]*dec(2,1ll*f[i]*g[i]%mod)%mod;
        ntt(g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) g[i] = 1ll*g[i]*Inv%mod;
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

inline void polySqrt(int* g,int* h,int n) {
    static int f[maxn],bit[30],b_len=0,B,lim,w,Inv,_g[maxn];
    while(n^1) bit[++b_len] = n, n = n+1>>1; g[0]=1;
    for(; b_len; --b_len) {
        w=bit[b_len], B=0, lim=1;
        while(lim<(w<<1)) lim<<=1, ++B;
        memcpy(f,h,sizeof(int)*w);
        memset(f+w,0,sizeof(int)*(lim-w));
        memset(_g,0,sizeof(int)*lim);
        polyInv(_g,g,w); // 没有预处理 rev[] 的原因是在 polyInv() 中已经处理过了
        ntt(f,lim), ntt(_g,lim);
        for(int i=0;i<lim;++i) _g[i] = 1ll*f[i]*_g[i]%mod;
        ntt(_g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) 
            g[i] = 1ll*iv2*(1ll*_g[i]*Inv%mod+g[i])%mod;
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

int main() {
    n=read(9);
    for(int i=0;i<n;++i) F[i]=read(9);
    polySqrt(G,F,n);
    for(int i=0;i<n;++i) print(G[i],' '); puts("");
    return 0;
}

3. 多项式求 \(\rm ln\)

3.1. 问题描述

给定多项式 \(f\)，求多项式 \(g\) 满足 \(g\equiv \ln(f)\ (\text{mod}\ x^n)\).

3.2. 解法

记 \(h(x)=\ln x\)，那么其实就是求

\[g\equiv h(f)\pmod{x^n} \]

两边同时求导得^[3]

\[g'\equiv \frac{f'}{f}\pmod{x^n} \]

在求出 \(g'\) 后，再求一次积分即可

\[\int x^a\ {\rm d}x=\frac{1}{a+1}\cdot x^{a+1} \]

另外还有一个需要注意的地方：

在模意义下 当且仅当 \([x^0]f(x)=1\) 时，\(f(x)\) 才有对数多项式，且常数项为零；当且仅当 \([x^0]f(x)=0\) 时，\(f(x)\) 才有指数多项式，且常数项为 \(1\).

解释一下为啥是这样（当然是感性的）。考虑 \(\ln f(x)\) 在 \(x_0=0\) 处进行泰勒展开时，第一项是 \(\dfrac{\ln f_0}{0!}\cdot x^0\)，若 \(f(x)\) 的常数项不为 \(1\)，首先 \(<1\) 肯定是未定义的，对于 \(>1\) 的情况，可以利用 \(\ln(x+1)\) 来观察，当 \(x>0\) 时这个函数在模意义下是不收敛的，所以没法算。\(\exp\) 也是同样的道理，这里就不讲了。

3.3. 代码实现

# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)

template <class T>
inline T read(const T sample) {
    T x=0; char s; bool f=0;
    while(!isdigit(s=getchar())) f|=(s=='-');
    for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
    return f? -x: x;
}
template <class T>
inline void write(T x) {
    static int writ[50], w_tp=0;
    if(x<0) putchar('-'), x=-x;
    do writ[++w_tp]=x-x/10*10, x/=10; while(x);
    while(putchar(writ[w_tp--]^48), w_tp);
}

# include <cstring>
# include <iostream>
using namespace std;

const int maxn = 4e5+5;
const int mod = 998244353, g = 3;

int inv(int x,int y=mod-2,int r=1) {
    for(; y; y>>=1, x=1ll*x*x%mod)
        if(y&1) r=1ll*r*x%mod; return r;
} const int ig = inv(g);
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }

int rev[maxn],F[maxn],G[maxn],n;

void ntt(int* f,int lim,bool opt=0) {
    int tmp, w, wn;
    for(int i=0;i<lim;++i) 
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1; mid<lim; mid<<=1) {
        wn = inv(opt? ig: g,(mod-1)/(mid<<1)), w=1;
        for(int i=0; i<lim; i+=(mid<<1), w=1)
            for(int j=0; j<mid; ++j, w=1ll*w*wn%mod) {
                tmp = 1ll*w*f[i|j|mid]%mod;
                f[i|j|mid] = dec(f[i|j],tmp),
                f[i|j] = inc(f[i|j],tmp);
            }
    }
}

void polyInv(int* g,const int* h,int n) {
    static int lim, B, bit[30], b_len=0, f[maxn], Inv, w;
    while(n^1) bit[++b_len]=n, n=n+1>>1; g[0]=inv(h[0]);
    for(; b_len; --b_len) {
        w=bit[b_len], lim=1, B=0;
        while(lim<(w<<1)) lim<<=1, ++B;
        for(int i=0;i<lim;++i)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<B-1);
        memcpy(f,h,sizeof(int)*w),
        memset(f+w,0,sizeof(int)*(lim-w));
        ntt(f,lim), ntt(g,lim);
        for(int i=0;i<lim;++i)
            g[i] = 1ll*g[i]*dec(2,1ll*f[i]*g[i]%mod)%mod;
        ntt(g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) g[i] = 1ll*g[i]*Inv%mod;
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

void dera(int* f) {
    for(int i=1;i<n;++i)
        f[i-1] = 1ll*f[i]*i%mod;
    f[n-1] = 0;
}

void inga(int* f) {
    for(int i=n-2;i>=0;--i)
        f[i+1] = 1ll*f[i]*inv(i+1)%mod;
    f[0]=0;
}

void polyLn(int* g,int* h,int n) {
    static int lim=1, Inv;
    while(lim<(n<<1)) lim<<=1;
    polyInv(g,h,n); dera(h);
    ntt(g,lim), ntt(h,lim);
    for(int i=0;i<lim;++i)
        g[i] = 1ll*g[i]*h[i]%mod;
    ntt(g,lim,1), Inv = inv(lim);
    for(int i=0;i<n;++i) g[i] = 1ll*g[i]*Inv%mod;
    inga(g);
}

int main() {
    n=read(9);
    for(int i=0;i<n;++i) F[i]=read(9);
    polyLn(G,F,n);
    for(int i=0;i<n;++i) print(G[i],' '); puts("");   
    return 0;
}

4. 多项式求 \(\rm exp\)

4.1. 问题描述

给定多项式 \(h(x)\)，求多项式 \(f(x)\) 满足 \(f(x)\equiv {\rm e}^{h(x)}\ (\text{mod}\ x^n)\).

4.2. 解法

由于之前已经解决了多项式求 \(\ln\)，所以不妨将这个玩意转化一下

\[\ln f(x)\equiv h(x)\pmod{x^n} \]

使用牛顿迭代来解决这个问题，记 \(g(f(x))=\ln f(x)-h(x)\)，由于它在模 \(x^n\) 意义下是为零的，所以有

\[\begin{align} f(x)&=f_0(x)-\frac{\ln f_0(x)-h(x)}{g'(f_0(x))}\\ &=f_0(x)-\frac{\ln f_0(x)-h(x)}{\frac{1}{f_0(x)}}\\ &=f_0(x)\cdot (1-\ln f_0(x)+h(x)) \end{align} \]

由于一切操作都是并行的，所以复杂度仍然是 \(\mathcal O(n\log n)\) 的！哦！让我们高声地赞美牛顿迭代法吧！~~常数大得飞起就是了~~。

4.3. 代码实现

# include <cstdio>
# include <cctype>
# define print(x,y) write(x), putchar(y)

template <class T>
inline T read(const T sample) {
    T x=0; char s; bool f=0;
    while(!isdigit(s=getchar())) f|=(s=='-');
    for(; isdigit(s); s=getchar()) x=(x<<1)+(x<<3)+(s^48);
    return f? -x: x;
}
template <class T>
inline void write(T x) {
    static int writ[50], w_tp=0;
    if(x<0) putchar('-'), x=-x;
    do writ[++w_tp]=x-x/10*10, x/=10; while(x);
    while(putchar(writ[w_tp--]^48), w_tp);
}

# include <cstring>
# include <iostream>
using namespace std;

const int maxn = 4e5+5;
const int mod = 998244353, g = 3;

int inv(int x,int y=mod-2,int r=1) {
    for(; y; y>>=1, x=1ll*x*x%mod)
        if(y&1) r=1ll*r*x%mod; return r;
} const int ig = inv(g);
int dec(int x,int y) { return x-y<0?x-y+mod:x-y; }
int inc(int x,int y) { return x+y>=mod?x+y-mod:x+y; }

int rev[maxn],F[maxn],G[maxn],n;

void ntt(int* f,int lim,bool opt=0) {
    int tmp, w, wn;
    for(int i=0;i<lim;++i) 
        if(i<rev[i]) swap(f[i],f[rev[i]]);
    for(int mid=1; mid<lim; mid<<=1) {
        wn = inv(opt? ig: g,(mod-1)/(mid<<1)), w=1;
        for(int i=0; i<lim; i+=(mid<<1), w=1)
            for(int j=0; j<mid; ++j, w=1ll*w*wn%mod) {
                tmp = 1ll*w*f[i|j|mid]%mod;
                f[i|j|mid] = dec(f[i|j],tmp),
                f[i|j] = inc(f[i|j],tmp);
            }
    }
}

void polyInv(int* g,const int* h,int n) {
    static int lim, B, bit[30], b_len=0, f[maxn], Inv, w;
    while(n^1) bit[++b_len]=n, n=n+1>>1; g[0]=inv(h[0]);
    for(; b_len; --b_len) {
        w=bit[b_len], lim=1, B=0;
        while(lim<(w<<1)) lim<<=1, ++B;
        for(int i=0;i<lim;++i)
            rev[i] = (rev[i>>1]>>1)|((i&1)<<B-1);
        memcpy(f,h,sizeof(int)*w),
        memset(f+w,0,sizeof(int)*(lim-w));
        ntt(f,lim), ntt(g,lim);
        for(int i=0;i<lim;++i)
            g[i] = 1ll*g[i]*dec(2,1ll*f[i]*g[i]%mod)%mod;
        ntt(g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) g[i] = 1ll*g[i]*Inv%mod;
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

void dera(int* f,int n) {
    for(int i=1;i<n;++i)
        f[i-1] = 1ll*f[i]*i%mod;
    f[n-1] = 0;
}

void inga(int* f,int n) {
    for(int i=n-2;i>=0;--i)
        f[i+1] = 1ll*f[i]*inv(i+1)%mod;
    f[0]=0;
}

void polyLn(int* g,const int* h,int n) {
    static int lim=1, Inv, f[maxn];
    while(lim<(n<<1)) lim<<=1;
    memset(g,0,sizeof(int)*lim);
    memcpy(f,h,sizeof(int)*lim);
    polyInv(g,h,n); dera(f,n);
    ntt(g,lim), ntt(f,lim);
    for(int i=0;i<lim;++i)
        g[i] = 1ll*g[i]*f[i]%mod;
    ntt(g,lim,1), Inv = inv(lim);
    for(int i=0;i<n;++i) g[i] = 1ll*g[i]*Inv%mod;
    inga(g,n); memset(g+n,0,sizeof(int)*(lim-n));
}

void polyExp(int* g,int* h,int n) {
    static int lim, B, bit[30], b_len=0, Inv, w, lng[maxn];
    while(n^1) bit[++b_len]=n, n=n+1>>1; g[0]=1;
    for(; b_len; --b_len) {
        w=bit[b_len], lim=1, B=0;
        while(lim<(w<<1)) lim<<=1, ++B;
        polyLn(lng,g,w); lng[0] = dec(lng[0],1);
        for(int i=0;i<w;++i) lng[i] = dec(h[i],lng[i]);
        ntt(lng,lim), ntt(g,lim);
        for(int i=0;i<lim;++i) g[i] = 1ll*g[i]*lng[i]%mod;
        ntt(g,lim,1); Inv = inv(lim);
        for(int i=0;i<w;++i) g[i] = 1ll*g[i]*Inv%mod;
        memset(g+w,0,sizeof(int)*(lim-w));
    }
}

int main() {
    n=read(9);
    for(int i=0;i<n;++i) F[i]=read(9);
    polyExp(G,F,n);
    for(int i=0;i<n;++i) print(G[i],' '); puts("");   
    return 0;
}

5. 多项式求幂

5.1. 问题描述

给定多项式 \(h(x)\)，求多项式 \(f(x)\) 满足 \(f(x)\equiv (h(x))^k\ (\text{mod}\ x^n)\).

5.2. 解法

先取对数

\[\ln f(x)\equiv k\ln h(x)\pmod{x^n} \]

再取一下指数

\[f(x)\equiv {\rm e}^{k\ln h(x)}\pmod {x^n} \]

5.3. 代码实现

Believe in yourself.

To be continued...

为了方便这里直接按无限展开来写，之后再说为什么。 ↩︎
这就是为啥可以直接无限展开，因为保留的项数是很少的。同时也可以发现我们对 \(g(x)\) 的项数存在一定限制。 ↩︎
这似乎是 \(\rm ln\) 的一个通用方法，将求 \(\rm ln\) 变成求逆元。 ↩︎

posted on 2022-05-28 15:22 Oxide 阅读(42) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部